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A  general  method  of  optimizing  continuous  microbial  growth  pro- 
cesses is  presented.  In  this  algorithm,  dynamic  information  is  used  to 
identify  parameters  of  a  dynamic  discrete  single-input/single  output 
model  equation.  The  steady-state  version  of  this  model  is  used  to 
predict  the  control  input  that  will  maximize  a  performance  index.  A 
major  advantage  of  this  method  is  no  requirement  of  a  priori  informa- 
tion. 

Numerical  simulations  on  a  Monod-type  model  are  shown.  In  these 
simulations,  biomass  productivity  was  maximized  with  respect  to  dilution 
rate,  temperature,  and  pH.  The  algorithm  successfully  found  the  optimum 
inputs  and  maintained  the  optimal  steady  states. 

An  anaerobic  digester  was  used  to  investigate  the  algorithm  on  an 
actual  process.  The  digester  used  was  a  6  liter  reactor  continuously 
fed  with  glucose  as  the  limiting  nutrient.  Digester  temperature^  was  the 
control  input  and  methane  gas  production  was  the  performance  index. 
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Two  innovations  were  introduced  to  improve  optimization.  The  first 
involved  the  use  of  an  alternate  sampling  basis  (as  opposed  to  sampling 
based  on  time)  to  speed  convergence.  In  this  method,  methane  production 
rate  was  measured  when  a  predetermined  amount  of  methane  was  collected 
(instead  of  time  elapsed).  The  second  innovation  used  a  constant  reac- 
tor yield  modus  operandi  (i.e.,  the  amount  of  substrate  fed  per  amount 
of  methane  produced  is  kept  constant)  to  reduce  the  chance  of  process 
failure  by  digester  imbalance.  These  changes  could  also  be  helpful  for 
other  microbial  growth  processes. 

Three  experimental  runs  were  made  on  the  digester,  two  in  the 
thermopnilic  range  and  one  in  the  mesophilic  range.  The  thermophilic 
runs  were  successful  at  increasing  gas  production  but  could  not  maintain 
optimal  production  because  of  the  sensitivity  of  the  methanogenic  popu- 
lation at  the  optimum  temperature.   The  mesophilic  run  optimized  at 
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CHAPTER  1 
INTRODUCTION 

In  the  past  score  of  years,  the  field  of  biotechnology  has  gained 
considerable  attention  by  providing  many  novel  processes  as  well  as 
alternatives  to  existing  chemical  processes.  Despite  unprecedented 
discoveries,  the  impact  of  this  field  in  industrial  applications  has 
been  severly  limited  due  to  the  complex  nature  of  using  living  organisms 
as  lilliputian  chemical  reactors.  Standard  procedures  used  in  the 
chemical  industry  simply  will  not  work.  The  development  of  new  and 
innovative  schemes  to  optimize  bioreactors  is  a  necessity  if  the  rich 
future  of  biotechnology  is  to  be  realized. 

In  the  past,  substantial  research  effort  has  been  expended  toward 
the  optimization  of  batch  and  semi-batch  bioprocesses  [1-7].  However, 
very  little  work  has  been  conducted  for  the  determination  of  the  steady- 
state  optimum  operating  conditions  of  continuous  bioreactors.  In  almost 
all  these  investigations  the  optimization  was  based  on  an  off-line 
identified  mathematical  model  or  on  ad  hoc  experimental  procedures.  The 
major  problem  of  all  methods  based  on  off-line  development  is  that  the 
process  conditions  vary  continually  due  to  population  shifts,  variable 
feed  concentrations,  wall  growth  and  other  disturbances.  Optimization 
based  on  these  models  is  questionable  in  most  cases. 

Consequently,  present  research  in  biotechnology  has  been  directed 
towards  the  on-line  determination  of  the  optimum  steady  state  of  contin- 
uous bioreactors  [8-12].  Traditional  solutions  to  the  on-line  optimiza- 
tion problem  involve  the  step  change  of  a  particular  manipulated 
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variable  and  the  evaluation  of  the  process  performance  measure  at  the 
subsequent  steady  state.  If  improvement  is  observed,  the  manipulated 
variable  is  changed  in  the  same  direction.  This  procedure  is  repeated 
until  no  further  improvement  in  the  performance  measure  can  be  ob- 
tained. The  obvious  drawback  of  such  methods  is  the  requirement  that  a 
steady  state  be  reached  at  each  step,  which  translates  into  a  long  time 
until  the  optimum  is  finally  achieved.  In  addition,  the  actual  process 
optimum  may  vary  at  a  faster  rate  than  it  can  be  located. 

To  enhance  the  speed  of  on-line  optimization  algorithms  for  contin- 
uous microbial  growth  processes,  dynamic  information  can  be  used  to 
estimate  parameters  of  a  dynamic  model.   The  optimal  constant  control 
value  is  determined  from  the  steady-state  version  of  the  identified 
dynamic  model.   Bamberger  and  Isermann  [13]  introduced  this  adaptive 
optimization  approach  in  1978.   In  their  work,  model  parameters  were 
identified  on-line  via  a  recursive  least  squares  estimator  and  a  Newton- 
Raphson  routine  was  employed  to  search  for  the  optimum.   This  approach 
has  been  used  on  many  processes  [9,11,12].   It  is  particularly  promising 
for  biochemical  processes  since  they  are  characterized  by  slow  dynamics. 
It  was  the  goal  of  this  work  to  develop  and  demonstrate  a  new 
method  for  the  adaptive  optimization  of  a  general  class  of  continuous 
bioreactors.   In  this  procedure,  a  model  with  parameters  identified  from 
on-line  data  is  used  to  explicitly  calculate  the  optimal  input.   Novel 
ideas  are  presented  in  the  areas  of  parameter  estimation  and  algorithm 
implementation.   In  addition,  a  constant  yield  operating  strategy  is 
proposed  to  improve  optimization  schemes. 

Chapter  2  gives  the  basic  background  and  approach  for  adaptive 
optimization   of   single-input/single-output   systems,   including   the 


problem  of  parameter  estimation.  The  specific  algorithm  is  presented  in 
the  following  chapter.  Numerical  simulations  are  presented  to  demon- 
strate the  ability  of  the  algorithm  to  successfully  locate  the  optimum 
temperature,  pH,  and  dilution  rate  for  biomass  productivity  in  continu- 
ous fermentors. 

In  the  following  chapters,  the  algorithm  is  applied  to  an  actual 
process,  anaerobic  digestion,  where  gas  production  is  optimized  with 
respect  to  temperature.  The  final  chapter  presents  the  conclusions  of 
this  work  and  offers  suggestions  for  future  work. 


CHAPTER  2 

GENERAL  ALGORITHM  DEVELOPMENT 

2.1.   Adaptive  Optimizafcion 
The  algorithm  presented  in  this  work  uses  data  acquired  on-line  to 
estimate  the  steady-state  value  of  a  control  variable  that  maximizes  a 
performance  index.   This  performance  index  can  generally  be  described  as 
a  function  of  the  measurement  output,  y,  and  the  control  input,  u: 


J=q(y(u) ,u) 


(2-1) 


Here,  y  is  indicated  as  a  function  of  u.  Optimum  steady-state  perform- 
ance  is  determined  when  a  constant  input,  u^,  satisfies  the  following 
conditions: 


dJ 


du 


dq(y(Ug),Ug) 


u  =u 
s  s 


du 


=  0 


u  =u 
s  s 


(2-2) 


d^J 


du 


<  0 


u  =u 

S   3 


(2-3) 


where  s  denotes  steady  state.   It  is  important  to  note  that  these  condi- 
tions do  not  specify  for  a  global  optimum.   Therefore,  if  the  system 

* 

exhibits  multiple  extrema,  the  input,  u^,  will  specify  only  for  a  local 


optimum,  which  may  or  may  not  be  global, 


In  order  to  evaluate  the  above  conditions,  it  is  necessary  to  have 
a  system  model  which  relates  the  output  to  the  input.  This  functiona- 
lity is  usually  written  in  continuous  form  as 

dy/dt  =  g(y,u;  6)  (2-4) 

where  9  is  a  vector  of  model  parameters  to  be  identified  on  line. 

However,  this  form  is  not  convenient  for  on-line  data  retrieval, 
where  measurements  are  taken  at  particular  times.  To  facilitate  com- 
puter implementation,  an  equally  valid  discrete-time  nonlinear  model  can 
be  employed. 

y(i  +  1)  =  f(y(i),  y(i-1),  ....  y(i-Ji),  u(i),  u(i-i),  ....  u(i-K);  G) 

(2-5) 

where  y(n)  and  u(n)  are  the  output  and  the  input  values,  respectively, 
at  the  nth  sampling  point. 

Model  structure  is  very  important  and  should  be  chosen  carefully. 
A  model  which  closely  approximates  the  actual  process  system  can  be 
employed  which  would  most  likely  entail  a  nonlinear  structure  with 
respect  to  parameters.  This  would  require  complex  search  methods  to 
estimate  9  and  the  predicted  optimum.  In  addition,  this  model  would 
require  a  priori  knowledge  of  the  process,  particularly  at  the  opti- 
mum. 

An  alternate  model  structure  can  be  chosen  to  simplify  parameter 
Identification  and  optimum  prediction.  The  most  convenient  model  would 
be  structured  in  such  a  way  as  to  estimate  a  unique  optimum  input,  i.e.. 


u*  should  have  only  one  value.   The  simplest  model  with  this  property 
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would  infer  the  steady-state  performance  index,  Jg,  has  a  parabolic 
shape  and,  therefore  is  quadratic  with  respect  to  Ug : 


2 

J  =  a  +  Yu  +  pu  (2-6) 

s         3     s 


where  a,    T,  and  p  are  functions  of  the  identified  parameters  6.   The 
evaluation  of  equation  (2-2)  yields: 


u*  =  -  Y/2p  (2-7) 
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with  the  following  sufficient  condition  from  equation  (2-3): 

p  <  0  (2-8) 

A  possible  drawback  to  this  line  of  reasoning  is  that  a  model  of 
the  type  described  above  would  in  general  be  only  locally  valid.  How- 
ever, this  does  not  present  a  major  problem  since  the  goal  of  the  adap- 
tive optimization  algorithm  is  not  necessarily  to  obtain  an  overall 
model  but  a  model  valid  at  the  locality  of  the  optimum.  Svoronos  and 
Lyberatos  [14]  proved  the  predicted  optimum  and  the  actual  process  op- 
timum would  coincide  if  the  model  structure  was  locally  valid  and  suffi- 
cient data  were  obtained  for  convergence  to  correct  parameter  esti- 
mates. Away  from  the  optimum,  the  algorithm  uses  the  locally  valid 
model  to  take  a  cautious  step  toward  the  predicted  optimum.  At  each  new 
location  new  data  can  be  acquired  to  update  the  parameters. 


For  example,  in  a  typical  case  where  the  measurement  output,  y,  is 
the  performance  index,  a  Hammerstein  type  model  [15]  can  be  employed  to 
meet  the  above  conditions: 

N  M  MM 

y(i)  =  I     a.y(i-j)  +  I     b.u(i-j)  +  Z   E  c  u(i-j )u( i-k)  +  d   (2-9) 
,j=1  ^  j  =  1  ^  j=1  l<  =  i  ^^ 

where  a.,  b.,  c.,  and  d  are  the  parameters,  and  N  and  M  are  the  dynamic 

J   J   J  *^ 

model  orders.  Before  such  a  model  can  be  implemented,  the  model  order 
and  sampling  period  must  be  appropriately  chosen.  Most  bioprocesses  are 
overdamped  systems  and,  therefore,  can  be  approximated  by  first  order 
systems  (N=1,M=1).  Higher  order  models  may  approximate  the  system 
closer,  but  they  will  also  increase  the  amount  of  data  needed  to  esti- 
mate parameters. 

The  choice  of  sampling  period,  T,  will  depend  upon  the  dominant 
time  constant  of  the  process.  If  the  value  of  T  is  too  large,  the 
measurement  output  at  each  interval  will  be  close  enough  to  the  steady- 
state  value  to  prevent  identification  of  a  dynamic  model.  The  T  chosen 
too  small  will  not  allow  the  measurement  output  to  change  significantly 
from  the  previous  measurement  and  could  lead  to  stability  problems  in 
parameter  estimation.  A  rough  estimate  of  T  is  usually  taken  as  ten  to 
twenty  percent  of  the  dominant  time  constant  [16].  However,  since  pro- 
cesses are  nonlinear,  the  dominant  time  constant  will  change  as  the 
process  moves  from  one  steady  state  to  another.  In  Chapter  5,  an  alter- 
nate sampling  basis  method  is  described  as  a  way  to  account  for  such 
problems. 


2.2.   Parameter  Estimation 

Ever  since  Gauss  and  Legendre  argued  over  who  first  developed  the 
method  of  least  squares,  a  veritable  plethora  of  parameter  estimators 
has  been  proposed.  However,  the  method  they  presented  in  the  early 
1800*3  still  remains  as  the  simplest  accurate  procedure  for  experimental 
systems  with  unknown  statistics.  In  i960,  Kalman  [17]  developed  a 
recursive  form  of  the  least-squares  method  and  supplied  a  practical 
procedure  for  digital  computer  implementation.  There  are  many  excellent 
texts  that  provide  the  basic  development  for  recursive  least  squares  as 
well  as  necessary  statistics  to  understand  the  problem  [18,19]. 

The  object  of  the  least  squares  is  to  minimize  the  square  of  the 
model  error.  The  error  for  measurement  i  is  defined  as 


e.  =  y(i)  -  hje  (2-10) 

1         -1- 


where,  referring  to  the  model  in  equation  (2-9)  as  an  example,  9  is  the 
vector  of  parameters 


^^y^2 ^N'^r^2 ''M'''ir''i2'''i3 "mm''^^ 


and  _h4  is  the  vector  of  previous  inputs  and  outputs 


[y(l-l),y(i-2) y(i-N),u(i-l ) ,u(i-2) , . . . ,u(i-N), 


u^(i-l ).u(i-1)u(i-2),u(i-1)u(i-3) u  (i-M),1], 


The  sum  of  the  squares  is  differentiated,  set  equal  to  zero,  and  solved 
for  e.  The  recursive  form  is  derived  by  considering  the  addition  of  new 
data  [20].  In  this  form,  new  parameter  estimates  can  be  calculated  with 
minimal  computations  and  less  storage  space  compared  to  batch  least 
squares.  A  problem  arises,  however,  with  recursive  least  squares  (RLS) 
when  large  amounts  of  data  are  acquired.  The  updating  becomes  insensi- 
tive to  new  data  so  that  the  new  parameter  estimates  are  largely  unaf- 
fected by  large  errors  in  the  model  equation.  To  overcome  this  problem, 
exponential  weighing  can  be  used.  This  can  be  expressed  in  recursive 
form  as: 

«..,(e)  =  A.,  _  1,   f9)  +  e^  (2-11) 


0  <  X  .  <  1  for  i  =  0  , 1  ,  .  . .  ,  iM- 1 

1 


where  2,,,  is  the  objective  function  and  X.    is  called  the  forgetting  fac- 
N  1 

tor.   The  updating  equations  are  obtained  by  differentiating  2...  with 

N 

respect  to  9  and  use  of  the  matrix  inversion  lemma  [18]. 


P.   h.  e. 

6,  =  9.  ,  +  ^"  "^   ' (2-12) 

1  +  h  P.  ,h. 
-1   i-1-i 

P.   h.h"!"  P. 
P.  =  7^  (P.  ,  -  '"^  :-"'   '"^)       (2-13) 

^    1        ^"^    1  +  h  P.  ,  h. 

-1   1-1  -1 


e^  =  y(i)  -  h[  6._^ 


(2-14) 


In  this  modification,  the  forgetting  factor  is  used  to  enlarge  the 
Govariance,  P,  thus  preventing  the  updating  of  the  parameters  from 
becoming  insensitive. 
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The  common  choice  of  \.  is  a  constant.  This  choice,  however,  has 
one  major  disadvantage.  As  parameters  converge  and  no  new  data  are 
acquired,  the  forgetting  factor  tends  to  "blow  up"  the  covariance  matrix 
[21].  To  correct  this  behavior,  Fortescue  et  al.  [22]  proposed  a  vari- 
able forgetting  factor  based  on  asymptotic  memory  length. 

2 

A.  =  1 =r-^ (2-15) 

(1  -  h,  P..^  h.)  a^ 

X.    =   \    .      if  X.  <  X  . 

1    min     1    min 

Here,  a  and  X  .   are  constants  dependent  upon  the  particular  system 
o      mm 

identified.  This  popular  method  has  come  to  be  known  as  RLSVFF  (recur- 
sive least  squares  with  variable  forgetting  factor).  Convergence  proofs 
and  other  characteristics  are  discussed  by  Cordero  and  Mayne  [23]. 

Recursive  least  squares  with  variable  forgetting  factor  is  used  as 
the  parameter  estimator  in  the  work  presented  in  the  following  sections 
with  two  innovations.  The  first,  described  in  Chapter  3,  involves  the 
inflation  of  the  covariance  matrix  for  sampling  periods  with  large 
errors.  In  Chapter  5,  a  method  is  presented  in  which  the  forgetting 
factor  is  applied  to  a  previous  covariance  matrix. 


CHAPTER  3 
NUMERICAL  SIMULATIONS 

3.1.   Inbroducblon 

Many  authors  have  proposed  algorithms  for  the  on-line  adaptive 
optimization  of  continuous  bioreactors  [8,10,12].  Whyatt  and  Petersen 
[11]  modified  the  previously  mentioned  approach  of  Bamberger  and  Iser- 
mann  [13]  to  include  optimal  input  trajectories  using  the  predicted 
optimum.  Simulations  on  a  contin'jous  fermentor  model  exhibited  success 
in  reducing  the  time  for  the  optimum  steady  state  to  be  reached. 

A  general  method  of  optimization  was  developed  by  Rolf  and  Lim 
[8,9]  using  a  gradient  scheme  to  find  the  optimum.  Their  method  was 
applied  to  maximize  the  biomass  productivity  of  a  continuous  culture  of 
baker's  yeast.  Experimental  results  were  obtained  by  Semones  and  Lim 
[24]  on  an  actual  fermenter  using  temperature  and  dilution  rate  as  the 
control  inputs.  Chang,  Pyjn ,  and  Lim  [25]  improved  the  parameter  esti- 
mation aspect  of  this  algorithm  by  introducing  a  bi level  forgetting 
factor.  Successful  experimental  and  simulation  studies  were  provided  by 
Chang  and  Lim  [26]. 

All  the  methods  listed  above  are  classic  examples  of  linear  adap- 
tive optimization,  i.e.,  a  linear  estimator  was  used  to  identify  param- 
eters of  a  linear  model  equation.  Golden,  Pangrle,  and  Ydstie  elected 
[12],  however,  to  address  the  problem  in  a  structured  format.  To  in- 
clude the  basic  static  nonlinearities  associated  with  continuous  biore- 
actors, they  proposed  the  use  of  a  nonlinear  least  squares  estimator  to 

adapt  a  structured  nonlinear  model  to  match  actual  operating  conditions. 

11 
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The  purpose  of  this  chapter  is  to  introduce  a  novel  adaptive  opti- 
mization scheme,  especially  suited  for  continuous  microbial  growth 
processes.  The  algorithm  implementation  is  easy,  and  the  use  of  gradi- 
ent schemes  is  not  required.  Furthermore,  no  previous  modeling  informa- 
tion is  required.  The  next  two  sections  describe  the  algorithm  in 
detail  and  provide  numerical  simulations  for  a  chemostat  model.  Biomass 
productivity  is  optimized  with  respect  to  dilution  rate,  temperature, 
and  pH. 

3.2.  Optimization  with  Respect  to  Dilution  Rate 
Since  the  major  goal  of  many  fermentation  processes  is  the  produc- 
tion of  biomass,  the  maximization  of  total  biomass  productivity  is  a 
major  design  objective.  A  suitable  performance  measure  should  include 
total  biomass  productivity  in  grams  of  cellular  material  per  unit 
time.  For  a  constant  volume  fermentation  the  product  Dx  where  D  is  the 
dilation  rate  and  x  is  the  biomass  concentration  is  clearly  such  a 
productivity  measure.  Furthermore,  since  the  major  running  cost  is  that 
of  the  utilized  substrate  a  reasonable  performance  index  is 

J  =  Dx  -  qDs°  (3-1  ) 

where  s"^  is  the  feed  substrate  concentration  and  q  is  factor  that  de- 
pends on  the  relative  cost  of  the  substrate  to  that  of  the  product, 
which  in  this  case  is  biomass. 

All  models  of  continuous  microbial  growth  processes  that  have 
appeared  in  the  literature  predict  that  the  above  performance  index  J(D) 
is  maximized  at  some  unique  value  of  the  dilution  rate  that  lies  between 
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zero  and  \i   ,    the  maximum  specific  growth  rate  at  the  prevailing  pH  and 
m 

temperature.  The  latter  two  operating  variables  for  the  p'irpose  of  this 
work  have  been  assumed  to  be  set  at  their  optimum  which  is  considered 
here  relatively  insensitive  to  process  changes.  The  existence  of  an 
optimum  dilution  rate  is  a  fact  that  has  also  been  repeatedly  verified 
experimentally  [27]. 

The  proposed  optimization  algorithm  determines  on-line  the  optimal 
dilution  rate.  In  addition,  it  ensures  that  the  fermentation  process  is 
operated  at  the  optimal  dilution  rate  at  all  times  thereafter,  despite 
changes  in  the  fermentation  environment. 

In  the  next  section  the  optimization  algorithm  is  described. 
Numerical  simulations  for  testing  the  algorithm  follow,  employing  con- 
tinuous chemostat  models  as  the  "true  processes." 

The  Adaptive  Optimization  Algorithm 

Since  fermentation  conditions  are  everchanging,  the  optimization 
algorithm  will  not  be  based  on  a  fixed  model;  instead  the  parameters  of 
the  model  will  be  continually  Identified  on  line.  But  before  parameters 
can  be  identified,  a  model  structure  must  be  assumed.  The  question  of 
what  structure  should  be  assumed  is  an  important  one. 

Fermentation  processes  are  nonlinear.  Furthermore,  the  type  of 
nonlinearities  involved  vary  depending  on  the  microorganisms  present  and 
the  fermentor  conditions,  and  as  a  result  are  not  generally  known.  It 
is  usual  to  assume  an  'unstructured  model  (e.g.  Monod  type)  whenever  a 
low  dimensional  model  is  needed,  but  these  models,  often  adequate  for 
steady-state  data,  fitting  generally  fail  to  match  the  dynamic  behavior 
of  the  system  [28].   Thus  it  would  be  dangerous  to  assume  a  particular 
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nonlinear  structure  and  base  a  general  adaptive  optimization  algorithm 

on  it. 

Nonlinear  systems  can  also  be  viewed  as  linear  systems  with  time- 
varying  parameters.  If  the  system  dynamics  are  slow,  these  parameters 
vary  slowly  in  time  and  as  a  result  can  be  identified  through  recursive 
parameter  estimators  which  forget  past  data  [22,29].  Fermentation 
processes  have  notoriously  slow  dynamics,  and  for  this  reason  a  meaning- 
ful linear  dynamic  model  can  be  continually  identifed  on  line,  regard- 
less of  what  the  "true"  nonlinear  model  of  the  process  is.  Thus  the 
optimization  algorithm  developed  here  is  based  on  a  linear  model. 

In  particular,  since  the  algorithm  is  to  be  implemented  using  a 
digital  computer,  the  system  is  modeled  via  the  linear  discrete-time 
input-output  form 


x(t+1)  =  a  x(t)  +  a,x(t-1)  +  •••  +  a  x(t-n) 
01  n 


+  b  D(t)  +  b,  D(t-1)  +  ...  +  b  D(t-m)  (3-2) 

o        1     .  m 


+  G 

Here  x(t)  denotes  biomass  concentration  at  time  x  (measured  in  n'omber  of 
sampling  intervals)  and  D(t)  denotes  dilution  rate  at  time  i  . 

The  above  dynamic  model  can  be  used  to  obtain  an  estimate  of  the 
optimal  steady-state  value  of  the  dilution  rate.  In  particular,  ne- 
glecting the  time-dependence  of  the  parameters,  model  (3-2)  gives  the 

following  pseudo  steady-state  equation: 

m 

I  ^ 

X  =   ^  =  Q  '         D  -  ^ (3-3) 

3       n      3        n 

1-  la.  1-  E  a 

i=0  ^  i=0 
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where  the  subscript  s  denotes  pseudo  steady  state.  Using  (3-3)  the  per- 
formance measure  (3-1 )  becomes  a  quadratic  in  D  and  can  be  solved 
explicity  for  the  maximizing  dilution  rate  D*  to  yield 


n 

c  -  qs  (1  -    a  ) 

D»(e)  =  i^U_J_  (3-4) 

m 

-2  Z  b. 
i=0  ' 


where  8^  is  a  vector  with  components  the  parameters  of  model  (3-2)  (see 
Appendix  B).  It  should  be  noted  that  only  equation  (2-2)  is  used  here 
to  designate  the  optimum.  The  optimallty  condition  of  equation  (2-3)  is 
not  included  because  it  was  not  needed  in  the  subsequent  numerical 
simulations.  This  was  mainly  due  to  the  lack  of  a  zero  order  term  in 
the  estimated  steady-state  performance  index  (i.e.,  J  =  0  when  D  = 
0).  Any  data  retrieved  with  J  >  0  and  D  >  0  would  tend  to  cause  the 
estimator  to  predict  a  properly  shaped  parabola.  In  addition,  for  this 
particular  Monod  type  system,  the  shape  of  x  plotted  versus  D  is  such 
that  good  parameter  estimates  would  predict  a  performance  index  which 
satisfies  equation  (2-3).  In  view  of  these  facts,  the  violation  of 
equation  (2-3)  is  highly  unlikely. 

However,  in  the  general  case  this  condition  of  optimality  should  be 
included  to  prevent  the  prediction  of  a  minimum  instead  of  the  expected 
maximum.  This  is  especially  true  under  the  following  three  circum- 
stances. First,  inflection  points  in  the  process  curve  (see  Figure  3-5) 
can  occur.  Second,  the  inclusion  of  a  zero  order  term  would  make  the 
predicted  performance  index  "free  floating,"  i.e.,  it  would  not  be  tied 
down  to  a  particular  point  (such  as  J  =  0  when  D  =  0).    Lastly, 
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measurement  noise  in  an  actual  system  could  temporarily  cause  poor 
parameter  estimates  resulting  in  poor  prediction. 

Since  the  parameter  vector  9^  is  not  known,  an  estimate  e(t)  is 
obtained  on  line  and  is  used  in  its  place.  Furthermore,  since  for  any 
fermentation  process  the  optimum  dilution  rate  cannot  exceed  1  hr  we 
limit  D*(e(t))  between  0  and  1.  The  basic  idea  of  the  algorithm  is  then 
to  change  at  every  sampling  instant  the  dilution  rate  D(t)  towards  the 
calculated  D  (£(t))  .  If  the  process  were  truly  linear  (i.e.  the 
parameters  of  model  (3-2)  constant),  this  algorithm  would  clearly  lead 
to  the  true  optimum  dilution  rate.  However,  since  the  true  process  is 
nonlinear,  the  success  of  the  above  scheme  is  not  evident.  In  a  pre- 
vious work  [30]  a  theorem  was  given  showing  this  to  be  indeed  the  case, 
no  matter  what  the  true  nonlinear  model  is,  provided  that  the  parameters 
of  the  local  linear  model  converge  to  the  correct  values. 

If  the  algorithm  sets  the  dilution  rate  to  the  predicted  optimum  at 
each  interval,  (D(t)  D*(a(t)),  large  differences  between  D(t)  and  D(t-1) 
may  appear.  The  true  parameters  would  as  a  result  vary  substantially 
with  time  leading  to  a  possible  failure  of  the  estimator  and  hence  of 
the  algorithm.  Furthermore,  since  the  original  parameter  estimates  will 
generally  be  poor,  taking  a  full  step  towards  the  calculated  optimum 
(D(t)  =  D*(9(t)))  could  lead  to  very  poor  performance  and  possibly 
washout.  These  possibilities  can  be  avoided  if  only  a  partial  step 
towards  D*(a(t))  is  taken,  i.e. 

D(t)  =  hD(t-l)  +  (1-h)  D»  (e(t));  0  <  h  <  1  (3-5) 
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This  equation  may  be  interpreted  as  first  order  filtering  of 
D*(0(t))  with  filter  parameter  h.  The  simplest  choice  for  h  is  a 
constant.  However,  in  order  to  ensure  safe  operation  when  the  parameter 
estimates  are  poor  (as  will  typically  be  the  case  at  the  beginning),  a 
conservative  choice  of  h  (i.e.  close  to  1)  would  have  to  be  made,  and 
this  would  result  in  a  very  slow  approach  to  the  optimum.  As  the  param- 
eter estimates  improve  with  time,  a  smaller  h  would  be  safe  and  would 
considerably  speed  up  the  convergence  of  the  algorithm. 

This  advocates  the  use  of  a  variable  h,  that  depends  on  the  quality 
of  the  parameter  estimates.  In  a  study  by  Harmon  et  al.  [30]  such  a 
variable  h  was  used  but  was  designed  in  a  rather  ad  hoc  manner. 

In  this  work  a  more  direct  and  meaningful  means  of  varying  h  is 
introduced.  Whenever  the  a  posteriori  estimation  error 


n     ^  ra  ^ 

v(t)  =  x(t)  -  Z  a,(t)  x(t-i-l)  -  I  b,(t)  D(t-i-l)  -  G(t)   (3-6) 
i=0  ^  i=0  " 


is  large,  the  parameter  estimates  are  poor  and  a  large  h  is  appropri- 
ate. Tnis  suggests  that  h  could  be  taken  to  be  a  raonotonically  increas- 
ing function  of  v(t).   A  reasonable  choice  is 


•^(t)  =    I  ^^^^1  (3-7) 

a   +  I  v(t)| 


where,  in  order  to  protect  against  a  low  v  due  to  cancellation  of  er- 
rors, we  filter  it  whenever  it  decreases,  i.e.  we  use: 


v(t)  if     v(t)  ^  v(t-l ) 

v(t)  =1  _  (3-8) 

(1-({.)v(t)  +  <pv{t-^)  if     v(t)  <  v(t-l) 
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Expression  (3-7)  contains  a  saturation  effect  similar  to  that  of  the 
Michael  is -Men  ten  rate  law  with  a  serving  the  role  of  a  saturation 
parameter.  The  lower  the  value  of  a  the  more  conservative  h  becomes. 
Equation  (3-7)  is  also  justified  by  the  arguments  given  in  Appendix  A, 
which  in  addition  provide  a  method  for  tuning  the  parameter  a  on  line. 

As  pointed  out  earlier,  the  success  of  the  proposed  algorithm 
depends  largely  on  the  ability  of  the  estimator  to  converge  to  the 
correct  parameter  values.   Since  the  biochemical  process  is  nonlinear 
and  the  model  is  linear,  this  requires  that  the  estimates  are  based  on 
local  data  only  (for  which  the  nonlinearities  are  negligible).   An 
estimation  algorithm  with  forgetting  discards  the  influence  of  past 
information  and  thus  the  estimates  are  mainly  based  on  recent  (and  hence 
local  for  a  slowly  varying  biochemical  system)  data.   However,  if  con- 
stant forgetting  is  employed,  the  estimation  algorithm  runs  into  prob- 
lems (e.g.  covariance  blowup)  whenever  it  stops  collecting  sufficient 
dynamic  information  (i.e.  in  the  neighborhood  of  a  steady  state).   The 
recursive  least  squares  estimator  with  variable   forgetting  factor 
(RLSVFF)  of  Ydstie  and  his  coworkers  [29,30]  solve  this  problem  by 
ceasing  to  forget  whenever  the  system  goes  to  a  steady  state,  a  desired 
effect  since  then  all  data  collected  are  local.   An  alternative  tech- 
nique is  to  periodically  reset  the  covariance  matrix.   It  was  found  that 
best  results  are  obtained  if  a  UDU  factorization  [31]  of  an  estimator 
that  is  essentially  a  combination  of  the  two  techniques  is  used  (see 
Appendix  B). 

A  problem  may  arise  if  the  step  size  h(t)  stays  close  to  1  for  long 
periods  of  time.  The  reason  is  that  D(t)  would  be  approximately  con- 
stant and  as  a  result  parameter  identifiabili ty  problems  could  arise. 
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To  avoid  this,  random  excitation  can  be  added  to  the  calculated  D(t). 
This  should  be  considerable  when  h(t)  is  high  but  very  low  when  the 
estimator  has  converged  (in  which  case  h(t)=0);  otherwise  it  would  then 
prevent  the  process  from  staying  very  close  to  the  optimum  steady 
state.  These  arguments  suggest  the  use  of  a  variable  excitation  level 
that  is  proportional  to  the  step  size  h(t)  of  equation  (3~7).  In  parti- 
cular we  set  the  excitation  at 

V(t)  =  (0.96  h(t)  +  0.04)  r(t)  (3-9) 

^OL        ''  '^'^    >  °TOL 

r(t)  =  {   w(t)        if       w(t)  e  C-D^q^'  ^TOL^ 

-^TOL        ''  "^'^    <  -^TOL 

Thus  V(t)  is  never  allowed  to  surpass  a  certain  tolerance  level  Di>g[_^, 
usually  chosen  about  20^  of  a  typical  dilution  rate  value.  The  white 
noise  w(t)  has  standard  deviation  D^q, /2,  which  implies  that  the  bound 
in  equation  (3-9)  is  active  less  than  5%  of  the  time.  Even  though  r(t) 
is  zero  mean  the  sum  of  r(t)  could  sometimes  become  significantly  dif- 
ferent from  zero  at  a  givf^n  time.  If  this  happens  when  h(t)  is  high 
then  we  are  not  exciting  the  process  about  the  desired  point.  To  avoid 
this  we  reverse  the  sign  of  r(t)  in  the  rare  case  when  the  absolute 
value  of  the  siam  of  r(t)  exceeds  2  D™q,  .  The  above  analysis  leads  to 
the  following  algorithm: 

STEP  1 :   Obtain  the  new  biomass  concentration  measurement 

x(t) 
STEP  2:   Use  the  estimator  of  Appendix  B  to  update  the 
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parameter  estimate  vector  e(t) 
STEP  3:   Calculate  the  step  size  via  equations  (3-6),  (A-9), 

(3-8)  and  (3-7). 
STEP  4:   Evaluate  the  desired  excitation  level  V(t)  using 

equation  (3-9). 
STEP  5:   Determine  the  estimated  optimum  dilution  rate 

0  ^  D»(e)  <  1  by  equation  (3-4). 
STEP  6:   Compute  and  implement 

D(t)  =  h(t)  D(t-1)  +  (1-h(t))D*  (9(t))  +  V(t) 
STEP  7:   Go  to  step  1 . 

Numerical  Simulations  for  Testing  the  Optimization  Algorithm 

In  order  to  test  the  behavior  of  the  developed  adaptive  optimiza- 
tion algorithm,  a  dynamic  model  of  the  chemostat  has  been  used.  Since  a 
typical  Monod-type  model  cannot  account  for  the  time-lag  observed  in  the 
dynamic  response  of  the  chemostat  [28]  a  time-delay  model  was  'osed 
instead.  The  model  introduced  earlier  [32]  is  described  by  the  follow- 
ing equations: 


M  zx 

m 

X  = Dx 

k  +  z 

3 


(3-10) 


1    P  sx 

^  =  -  Y  •  ;r^  *  °^^  -^^ 

s 


z  =  a(s-z) 


where   \i     is    the  maximum  specific   growth  rate,    k^    is    the  Monod   constant, 
m  "^  °  '      s 

Y     is     the    yield     and     z     is     a    weighted    average    of    previous    substrate 
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concentrations.  The  parameter  a,  called  the  adaptability  parameter,  is 
a  measure  of  the  organism's  ability  to  adjust  its  growth  rate  when  a 
change  in  the  conditions  of  the  chemostat  is  brought  about. 

The  model  parameters  used  for  the  simulations  are  tabulated  in 
Table  3-1.  Parameters  for  the  optimization  algorithm  are  given  in  Table 
3-2.  Even  though  the  true  model  order  is  three,  a  second  order  model  of 
the  form  of  equation  (3~2)  is  used  in  order  to  test  the  robustness  of 
the  algorithm.  The  sampling/control  interval  is  five  minutes.  Initi- 
ally the  chemostat  is  operated  for  1 .75  hours  at  the  nonoptimum  dilution 
race  of  0.05  hrs"  ,  during  which  period  estimation  but  no  control  oc- 
curs. After  this  initializatioa  period,  the  full  optimization  algorithm 
is  Implemented.  After  100  hours  a  load  change  in  the  feed  substrate 
concentration  is  introduced  as  s^^  changes  from  30g/l  to  20g/l.  This 
shifts  the  optimum  value  of  the  performance  measure  from  J  =  2.22  to  J 
=  1.12.  A  second  disturbance  is  introduced  at  200  hrs;  a  change  in  tem- 
perature or  pH  has  been  assumed  to  cause  a  change  in  u  from  0.7  to  0.9 
hr~  .   This  causes  a  further  change  in  J  from  1.12  to  1.44. 

The  behavior  of  the  algorithm  is  shown  in  Figures  3-I  and  3-2. 
Figure  3-1(a)  shows  the  implemented  control  D(t).  The  calculated  opti- 
mum dilution  rate,  D*(6(t)),  is  shown  in  Figure  3-1 (b)  and  the  stepsize 
h  is  shown  in  Figure  3-1(c).  From  this  f"ig'are  it  is  obvious  that  when  a 
disturbance  is  introduced,  h  increases  rapidly  rendering  the  control 
behavior  more  conservative.  As  the  parameter  estimates  improve,  h  drops 
towards  zero,  and  the  excitation  level  progressively  decreases. 

Figure  3-2  shows  the  behavior  of  the  biomass  concentration  (Figiare 
3-2(a))  and  the  performance  measure  J  as  a  function  of  time  (Figure  3- 
2(b)).   On  the  same  fig'jre  a  dashed  line  shows  the  actual  optimum  values 
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Table   3.1.      Model  parameters    (at   time  t=0) 

u  =   0.7  hr~'^  S  =   30  g/X,  D   =  0.05  hr"'' 

m  o 


Y   =  0.5  a  =   3  hr"''  x  =    14.153  g/«.t 


k  =   22  g/ilt  q  =   0  3   =   z  =    1.6923   g/2.t 

3 


Table  3.2.  Algorithm  parameters 

T  =  5  min     (sampling  interval) 

n  =  m=2       (order  of  model  in  equation  2) 

(p   =  0.9       (filter  parameter  in  equation  8) 

D^QL  =  0.05  hr"^  (excitation  tolerance  level  (eq  9)) 

5  =  0.05      (limit  to  the  erroneous  part  in  D(t)  appearing  in 

equation  A4) 
0=10       (estimator  parameter) 

0 

T  =  10        (estimator  parameter) 
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"igure  3-1.  Temporal  variation  of  (a)  implemented  dilation  rate,  (b) 
calculated  opti.TiUm  cilition  rate,  both  in  (hours),  and  (c) 
step  size  using  time  delay  model  as  true  process.  Distur- 
bances occur  at  '00  and  200  hours. 
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of  the  performance  measure.  It  is  seen  that  the  adaptive  optimization 
algorithm  successfully  locates  the  optimum,  even  after  serious  dist'or- 
bances  in  the  chemostat  environment  are  introduced  changing  its 
position. 

In  order  to  establish  the  generality  of  the  proposed  algorithm,  its 
success  was  established  with  different  models  for  continuous  bioreac- 
tors,  including  models  with  substrate  inhibition  and  variable  yield 
factor.  The  algorithm  was  successful  in  taking  all  these  processes  to 
the  optimum  without  having  to  make  any  changes  in  the  tuning  parameters. 

Of  particular  interest  is  the  behavior  of  the  algorithm  for  a 
process  model  taken  from  the  literature  [33].  This  r-ommercially  impor- 
tant process  is  the  production  of  single  cell  protein  using  methylo- 
trophs.  The  successful  behavior  of  the  algorithm  is  portrayed  in  Fig- 
ures 3-3  and  3-''^. 

3.3.  Optimization  with  Respect  to  Temperature  and  pH 
The  purpose  of  this  section  is  to  present  an  extension  of  the 
previous  section  that  can  be  applied  to  the  determination  of  optimal 
temperature  and  pH  values  for  a  continuous  fermentor  operating  at  steady 
state.  One  important  alteration  to  the  previous  dilution  rate  algorithm 
involves  the  form  of  the  identified  model.  Since  the  algorithm  requires 
an  explicit  formula  for  the  optimal  control  value  (something  not  pos- 
sible with  a  simple  linear  model),  a  nonlinear  model  of  the  Hammerstein 
form  is  used.   An  additional  modification,  which  leads  to  improvement  in 
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Figure  3-2.  Temporal  variation  of  (a)  biotnass  concentration  and  (b) 
performance  inaex  (dashed  line  represents  actual  optimum 
steady  state).   Process  model  as  in  Figure  3-1. 
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•igure  3-3.  Temporal  variation  of  (a)  implemented  dilution  rate,  (b) 
calculated  optimum  dilution  rate,  both  in  (hours),  and  (c) 
step  size  using  variable  yield  methylotroph  model. 
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Figure  3-'^.  Temporal  variation  of  (a)  biomass  concentration  and  (b) 
performance  index  (dashed  line  represents  actual  optimum 
steaay  state).   Process  model  as  in  Figure  3-3. 
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the  previous  algorithm  but  becomes  mandatory  for  the  case  of  temperature 
optimization,  is  the  monitoring  of  the  second  derivative  of  the  perform- 
ance measure  calculated  from  the  identified  model.  This  ensures  that 
the  calculated  control  value  (i.e.,  point  of  zero  slope)  results  in  a 
maximum  index  of  performance  instead  of  a  minimum.  The  calculation  of  a 
minimiam  occurs  sometimes  when  the  parameters  are  estimated  in  regions 
away  ^rom  the  actual  optimum. 

In  the  next  section,  the  algorithm  is  described.  Examples  of 
numerical  simulation  results  employing  a  Monod  model  with  time  delay  as 
a  true  process  are  presented  in  the  following  section. 

The  Adaptive  Control  Algorithm 

The  algorithm  presented  in  this  section  adaptively  estimates  the 
value  of  the  control  variable,  either  temperature  or  pH,  that  maximizes 
a  performance  index.  A  typical  performance  objective  is  the  production 
of  biomass  in  the  steady-state  operation  of  a  continuous  bioreactor. 
Such  a  measure  would  include  a  term  DXg  where  D  is  the  dilution  rate  and 
Xg  is  the  steady-state  biomass  concentration,  to  accoiont  for  producti- 
vity. Also,  since  a  major  operating  cost  is  for  the  utilized  substrate, 
a  reasonable  performance  index  is 


J  =  DX   -  qDS  (3-11  ) 

s      o 


where  S^  is  the  inlet  substrate  concentration  and  q  is  a  constant  de- 
pendent on  the  relative  cost  of  substrate  to  biomass.   However,  since  X 
is  the  only  terra  in  (3-11)  dependent  on  temperature  and  pH,  the  optimi- 
zation of  biomass  concentration  will  result  in  the  optimization  of  J. 
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To  obtain  the  optimum  control,  a  function  relating  X„  to  the  control 
variable  U„  is  needed. 

O 

In   many   cases   a   discrete-time   linear   model   can  be   used    to  approxi- 
mately relate   the  control   LI(t)    to   the   measured  output  X(t): 


X(t+1)   =  a^  X(t)    +  a^    X(t-1)   +    +   a^  X(t-n) 

+  bg  U(t)   +   b^    U(t-l)   +    +  b^  U(t-ra)   +  d  (3-12) 

Here,  t  denotes  time  in  number  of  sampling  intervals.  The  above  equa- 
tion relates  biomass  concentration  at  time  t+1  to  previous  biomass 
concentrations  and  previous  values  of  the  control.  This  model  was 
successfully  used  in  previous  work  [10]  to  determine  the  optimum  of  J 
with  respect  to  dilution  rate.  However,  if  U  is  temperature  or  pH  the 
performance  measure  is  optimized  at  a  bound  and  not  at  a  point  of  zero 
slope  (dJ/dUg  =  0).  Since  the  true  fermentation  optimum  is  at  a  point 
of  zero  slope,  a  model  with  this  feature  would  be  desirable.  In  related 
work  [14],  it  was  proven  that  if  such  a  model  has  the  correct  local 
parameters  then  the  model  optimum  is  identical  to  the  true  process 
optimum.  A  simple  way  of  constructing  such  a  model  is  to  modify  equa- 
tion (3-12)  to  a  Hammerstein  form  by  adding  terms  quadratic  in  U: 


X(t+1)  =  agX(t)  +  +  a^X(t-n)  + 


bQU(t)  +  +  bjjjU(t-m)  + 


c^^u2(t)  -  ....  *  c^u2(t-m)  ^  d  (3-13) 
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It  should  be  noted  that  the  large  majority  of  fermentors  are  overdamped 
and  therefore  a  first  order  model  (n=m=0)  could  be  adequate. 

Optimization  now  leads  to  an  explicit  analytical  solution  for  the 
optimal  control  value.  In  particular,  the  steady-state  version  of  (3- 
13),  (i.e.,  when  X(t)  =  X(t+1)  and  U(t)  =  U(t+1)  for  all  t),  yields  the 
control  value  for  which  J  is  maximized  (i.e.,  dJ/dUg  =  0): 


"opt  =  -^/2c  (3-14) 


m  mm 

where  b  =  E  b .   and   c  =  E   L     c .  . 

i  =  0  ^  i=0  j=i  ^^ 


Since  the  true  parameters  of  model  (3-15)  are  not  known,  they  are  esti' 
mated  from  sampled  data.   The  estimated  parameters,  denoted  by 


9=[a,a,,....,a,b,....,b,c   ,...., c  ,d] 
0I'     'no'     'moo'     'mm' 


are  updated  after  each  sampling  interval  via  a  recursive  least  squares 
routine  similar  to  the  algorithm  proposed  by  Fortescue,  Kershenbaum  and 
Ydstie  [22].  This  routine  uses  a  variable  forgetting  factor  to  discard 
old  data  which  may  bias  the  estimates,  6.  This  is  done  to  ensure  that 
the  parameters  are  estimated  from  local  data.  The  discarding  of  data 
ceases  when  the  system  converges,  since  all  acquired  data  are  then 
local.  In  addition  to  forgetting,  a  mechanism  for  covariance  matrix 
resetting  was  included.  This  allows  the  estimator  to  discard  a  larger 
number  of  data  when  disturbances  occur.  The  complete  set  of  estimator 
equations  can  be  found  in  Appendix  B. 
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The  parameter  vector,  6,  estimated  from  the  least  squares  routine 
is  used  to  calculate  the  estimated  optimal  control  value,  U^  l^,  from 
equation  (3~1^).   In  addition,  if  bounds  are  known  they  are  included,  so 

that 

U       if  -b/2G  >  U 
max  max 

U  ^  =   -b/2c    if  U  .   <  -b/2c  <  U  (3-15) 

opt  rain  max 

U  .      if  -b/2c  <  U  . 
rain  rain 

It  must  be  stressed  that  the  validity  of  the  calculated  control  depends 

directly  on  the  validity  of  the  estimated  parameters.   Therefore  the 

control  value  should  be  implemented  with  some  measure  of  caution.   In 

other  words,  a  step  toward  U  ,  should  be  taken  at  each  interval  instead 

opt 

of  implementing  U  =  U 

opt 

To  achieve  this,  a  first  order  filter  can  be  used.  The  control  law 
then  becomes: 


U(t)  =  (l-h)U  ^  +  hU(t-l);   0  <  h  <  1  (3-16) 

opt 


where  h  is  the  filter  parameter.  The  simplest  choice  for  h  is  a  con- 
stant. This  constant  would  have  to  be  chosen  conservatively  (close  to 
1 )  to  avoid  large  steps  in  the  wrong  direction  when  parameter  estimates 
are  poor.  However,  as  the  parameters  improve,  employing  a  large  h  will 
unnecessarily  slow  down  convergence  of  the  algorithm.  Therefore,  a 
variable  filter  parameter  based  on  estimation  errors  is  advocated.  As 
in  the  previous  section,  h  is  chosen  as  a  monotonically  increasing 
function  of  v(t),  the  one-way-filtered  a  posteriori  output  prediction 
error : 
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h(t)   -        l^^'^l  (3-17) 

a  +    |v(t)| 


v(t)  if  v(t)   ^   v(t) 

v(t)     =  (3-18) 


(1-$)v(t)   +  (t.v(t-1)   if  v(t)   <   v(t) 


n 
v(t)   =  X(t)   -     1     a,(t)   X(t-l-i) 


i=0 


m     ^ 

Z     b.(t)  U(t-l-i)  (3-19) 

i=0     ' 


ra       m       ^  P 

lie. .(t)    U   (t-1-i)   -   d(t) 
i=0  j=i       'J 


The  a  posteriori  error  v(t)  is  filtered  whenever  it  decreases  in  order 
to  protect  against  low  values  due  to  cancellation  of  errors  in  the 
parameter  estimates.  The  parameter  a  is  analogous  to  a  saturation 
parameter  in  a  Michaelis-Menten  type  rate  expression.  The  lower  is  its 
value,  the  more  conservative  is  the  filter  in  (3-16).  Using  arg<jments 
presented  in  section  3.2,  the  following  formula  for  t'oning  a  on  line  can 
be  derived 

^  ^      26^U^(t-1)|c(t)|  (3_20) 

|2  U      ,  /U(t-1  )    -^    1  I 
'        opt  .  ' 

The  parameter  6  can  be  thought  of  as  an  upper  bound  to  the  relative 
error  in  U(t).  For  example,  if  U  is  temperature  of  value  in  the  order 
of  300  K  a  choice  6  =  .0001  limits  the  erroneous  part  in  control  law  (3- 
16)  to  0.03  K. 
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A  white  noise  term  is  added  to  the  right  hand  side  of  equation  (3- 
16)  for  excitation.  It  is  included  to  prevent  U(t)  from  remaining 
constant  for  long  periods  of  time  thus  causing  parameter  identifiability 
problems.  This  terra  is  designed  for  maximum  excitation  at  h=1  (i.e., 
poor  parameter  estimates).  However,  when  the  estimator  has  converged 
(in  which  case  h  is  very  low),  excitation  should  be  low  to  avoid  keeping 
the  process  far  from  the  optimum  steady  state.  Consequently,  the  imple- 
mented excitation  is  taken  as: 

V(t)  =  (0.96  h(t)  +  O.Oi|)*r(t)  (3-21) 

Utol  when  w(t)  >  Utol 
r(t)  =    w(t)  when  w(t)  e  [-Utol,  Utol]        (3-22) 
-Utol  when  w(t)  <  -Utol 

where  w(t)  is  white  noise  with  zero  mean  and  variance  equal  to  Utol/2. 
Thus  the  absolute  value  of  V(t)  is  never  allowed  to  surpass  the  toler- 
ance level,  Utol,  usually  chosen  to  be  about  ^0%  of  the  control  variable 
range.  Tolerance  levels  are  employed  to  prevent  large  movements  in  the 
control  variable.  Furthermore,  the  sign  of  r(t)  is  reversed  whenever 
the  absolute  value  of  the  s'Jin  of  V(t)  over  all  time  exceeds  2  Utol. 
This  is  done  to  ensure  that  the  control  variable  does  not  travel  over 
long  periods  of  time  when  h  is  high. 

To  complete  the  algorithm  one  more  feature  is  needed.  Since  bio- 
mass  concentration  should  not  be  minimized  with  respect  to  temperature 
or  pH  at  a  point  of  zero  slope,  it  is  expected  that  points  of  zero  slope 
are  maxima  and  therefore  the  estimated  second  derivative  satisfies 
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3U  ' 


2g 

1  -  a 


<  0 


(3-23) 


However,  this  might  not  be  true  due  to  convex  o'jrve  shapes  or  to  estima- 
tion errors.  In  this  case  equation  (3-1 '^)  should  not  be  used  since  it 
would  lead  towards  a  predicted  minimum.  Instead,  the  control  variable 
should  be  moved  in  the  opposite  direction.  Since  the  violation  of  (3- 
23)  may  be  due  to  inaccurate  parameter  estimates  the  control  step  should 
be  small  and  excitation  V(t)  should  be  added  to  aid  parameter  estima- 
tion. However,  the  step  size  should  not  be  taken  so  small  as  to  allow 
the  excitation  to  reverse  the  direction  of  control  movement.  Thus  it  is 
taken  equal  to  the  maximum  excitation  level  Uf-^j^  (see  equation  (3- 
22)).  Of  course  if  the  parameter  estimates  are  known  to  be  worthless 
(h  =  1  )  no  trust  should  be  placed  on  the  sign  of  the  estimated  second 
derivative  of  J.  In  this  case  the  control  should  only  be  excited  lantil 
the  parameter  estimates  improve  (e.g.  h  drops  below  0.99). 

In  summary  the  above  analysis  leads  to  the  following  algorithm: 


STEP  1 : 
STEP  2: 

STEP  3: 

STEP  4: 
STEP  5A: 


Obtain  the  new  biomass  concentration  measurement 

Use  the  estimator  of  section  3-2  to  update  the  parameter 

estimates  9 

Calculate  the  filter  parameter  h(t)  via  equations  (3-17)-(3- 

20) 

Check  i^  inequality  (3-23)  is  satisfied 

If  inequality  (3-23)  is  satisfied  calculate 

U  ^  via  equation  (3-15) 
opt 

V(t)  via  equations  (3-22)  and  (3-21) 
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and  implement 

U(t)  =  h(t)U(t-1)  +  [1-h(t)]  U  ^  +  V(t) 

opt 

STEP  5B:   If  inequality  (3-23)  is  violated  and  h  <  0.99  take  a  step  away 

from  U   ^ ,  i.e. 
opt 

U(t)  =  U(t-1)  -  sgnL-b/2c-U(t-1 )]  U^  ,  +  V(t) 

tol 

where  V(t)  is  obtained  from  equations  (3-22)  and  (3-21) 
STEP  5C:  1^     inequality  (3-23)  is  violated  and  h  >  0.99  excite  the 

control  variable,  i.e. 

U(t)  =  U(t-1  )  *   V(t) 
STEP  6:    Go  to  step  1  . 

The  above  described  algorithm  has  several  parameters  that  must  be 
set  by  the  user.  Some  trial  and  error  may  be  needed  and  it  is  hoped 
that  the  following  discussion  will  aid  the  task. 

The  sampling  interval  should  not  be  larger  than  10?  of  the  dominant 
time  constant  of  the  fermentor.  There  is  an  advantage  to  going  to  lower 
values  since  it  would  tend  to  speed  up  convergence  of  the  algorithm. 
However,  the  sampling  period  value  should  never  be  lowered  to  the  point 
that  under  dynamic  conditions  successive  measurements  do  not  differ 
within  the  precision  of  the  measuring  device. 

The  estimator  used  has  two  parameters.  One  of  them,  T,  does  not 
significantly  affect  the  algorithm  performance  as  long  as  it  is  set  to  a 
low  enough  value  (less  than  10"  ).  The  other  parameter,  o  ,  adjusts  the 
rate  at  which  the  estimator  forgets  past  information.  The  higher  o  is, 
the  less  is  the  forgetting;  and  the  less  the  forgetting,  the  more  the 
steady-state  offset  that  can  result  from  the  algorithm  [T4].  On  the 
other  hand,  excessive  forgetting  from  a  very  low  value  of  o  could  cause 
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estimator  failure.   This  would  be  reflected  in  h(t),  which  would  contin- 
ually stay  very  close  to  1.   One  can  set  o  as  follows.   Start  from  a 

o 

very  conservative  value,  say  lO-'r,  where  r  is  the  expected  measurement 

variance  (see  [22]),  and  then  decrease  it  until  the  estimator  starts 

having  difficulties  (observe  h(t)).   In  most  cases  the  estimator  will 

perform  well  for  a  wide  range  of  a  values  (2-3  orders  of  magnitude). 

o 

The  parameter  (f  in  equation  (3-19)  affects  how  fast  h  drops  (the 
lower  it  is,  the  faster  h  decreases).  A  value  in  the  range  .7  -  .95  is 
recommended.  A  value  <t>  =  .9  implies  that  h  decreases  by  ^0%  (  l-<j)  )  as 
a  result  of  negligible  estimation  errors  at  5  consecutive  sampling 
instants.  Clearly,  the  larger  the  sampling  period  is,  the  lower  must  0 
be  to  achieve  the  same  rate  of  decrease  in  h.  Another  parameter  that 
affects  h  is  5  (equation  (3-20));  che  higher  it  is,  the  faster  h  de- 
reases.  Thus  the  larger  the  sampling  period  is,  the  higher  must  6  be 
to  obtain  the  same  rate  of  decrease  in  h.  As  mentioned  previously,  one 
can  be  guided  towards  choosing  6  by  thinking  of  it  as  a  measure  of  the 
relative  error   in  U(t). 

Finally,  the  parameter  U^^^  of  equation  (3-22)  does  not  depend  on 
the  sampling  period  and  the  previously  mentioned  value  {}0%  of  the 
control  variable  range)  should  be  satisfactory. 

Numerical  Simulations  for  Testing  the  Optimization  Algorithm 

In  this  section,  numerical  simulations  are  presented  to  test  the 
optimization  algorithm.  To  account  for  dynamic  behavior  in  a  chemostat, 
a  time-delay  Monod  model  was  employed  as  the  "true  process".  This  model 
was  introduced  by  O'Neil  and  Lyberatos  [32]  and  is  described  by  the 
following  equations: 
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M  ZX 
3 


,  P  sx 

1  — !5 —  +  D(s°-s) 

Y  K   +  S         ^ 
s 


Z  =  a(S-Z)  (3-2M) 

where  K  is  the  Monod  constant,  Y  is  the  yield,  X  is  biomass  concentra- 
tion, S  is  substrate  concentration  and  Z  is  a  weighted  average  of  previ- 
ous substrate  concentrations^  The  delay  parameter  a,  called  the  adapta- 
bility parameter,  is  a  measure  of  the  microbe's  ability  to  adjust  its 

growth  rate  to  changing  substrate  concentrations.   In  addition,  n   ,    the 

m 

maximum  specific  growth  rate,  is  taken  to  be  a  function  of  pH  and  tem- 
perature : 


^m   ^max  pH  T 

3,  T  e-^^»^ 

^T  "  ~  AS~7r   -AH,/RT  (3-25) 

1  +  e   d   e   d 

f  „  =  ^ (3-26) 

P"   (1  +  K^/LH*]  +  Ch']/k^) 

where  f  ^^  and  fj   are  enzyme  activity  functions  [3^^].  The  constants 

8   and  B  are  proportionality  factors,  while  Ki  and  Kp  are  equilibrium 

pH        1  '        '- 

constants  for  enzyme  ionization.  The  activation  energy  for  an  enzyme 
catalyzed  reaction  is  denoted  by  E,  whereas  AS  and  AH  denote  the 
entropy  and  enthalpy  of  deactivation,  respectively.  Values  for  the 
parameters  in  equations  (3-24)  -  (3-26)  are  given  in  Table  3-3. 
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Table  3-3.   Simulation  parameters 


Model 


y    =  0.7  hr'"^  3„  =  30  g/S,  D  =  0.05  hr~^ 

max  o 

Y    =  0.5  a  =  3  hr"''  K   =  22  g/«.lt-hr 


pH  Run  Initialization 


pH   =5.5         X  =  13.29         S  =  Z  =  3.'+1  T  =  312.0  K 


Temperature   Initialization 


Tq      =300  X      =    12.67  S  =   Z   =    J4.66  pH   =  7.0 


Enzyme   Activity 


gmol  d  gmole-K  d  gmol 


X^      =    10"^-^  K2     =    10"^"^ 


5^     =  9.845  X  lo"^  6  ^  =    I.O63 

1  pH 
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The  dependence  of  steady  state  biomass  concentration  on  temperature 
and  pH  is  shown  in  Figure  3-5.  The  optima  are  7.0  and  312.0  K  for  pH 
and  temperature,  respectively. 

The  numerical  simulations  were  performed  with  a  sampling  period  of 
five  minutes.  The  dilution  rate  of  the  chemostat  was  set  at  0.05  hr~'^ 
for  both  temperature  and  pH  optimization.  Initially,  the  control  vari- 
able was  set  at  a  nonoptimum  value  and  excited  with  white  noise  r(t) 
around  that  point  for  a  period  of  1.75  hours.  After  this  initialization 
period,  the  algorithm  was  implemented  in  full  with  parameters  as  given 
in  Table  3-4. 

The  optimization  with  pH  as  the  control  variable  is  given  in  Figiore 
3-6.  Initially,  the  pH  was  5.5  and  it  converged  to  the  optimum  with  an 
error  of  0.7%.  The  calculated  optimum,  equation  (3-14),  converged  to 
the  correct  value  within  4  hours  but  was  not  fully  employed  until  the 
filter  parameter,  h,  decreased. 

The  temperature  simulation,  shown  in  Figure  3-7,  was  initiated  at 
300  K.  The  temperature  subsequently  converged  to  within  0.1$  of  the 
true  optimum.  Initially,  the  estimator  identi*'ies  parameters  that 
result  in  a  minimum  instead  of  a  maximum  of  J.  Condition  (3-23),  there- 
fore, was  violated  and  Step  5B  was  implemented  to  drive  the  system  away 
from  this. 

In  conclusion,  the  adaptive  optimization  has  proven  very  successful 
in  bringing  a  simulated  fermentation  process  to  its  optimum  conditions 
and  ensuring  optimal  operation  thereafter.  The  algorithm  should  prove 
very  useful  in  two  types  of  application:  (i)  the  fast  determination  of 
optimal  temperature  and  pH  for  microbial  cultures  and  (il)  securing 
optimal  operation  for  large  scale  fermentation  processes. 
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Figure  3-5.   The  dependence  of  steady-state  biomass  concentration  on  (a) 
pH  and  (b)  temperature  using  equations  (3-2^4)  to  (3-26). 
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Fig-ure  3-6,  pH  optimization.  Temporal  variation  of  (a)  the  actual  pH 
control  value,  (b)  the  parameter  estimates,  and  (c)  step 
size. 
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Figure  3~7.  Temperature  optimization.  Temporal  variation  of  (a)  the 
actual  temperature,  (b)  the  predicted  optimum,  both  in  K, 
and  (c)  step  size. 
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Table  3-^^.   Algorithm  parameters 


(sampling  interval) 

(order  parameters  in  equation  3-13) 

(filter  parameter  in  equation  3-1 8) 

(parameter  of  equation  (3-22)  for  temperature 

simulation) 
(pH  simulation) 
(parameter  of  equation  (3-20)  for  temperature 

(pH  simulation) 

(forgetting  factor  parameter) 

(initial  covariance  matrix  parameter) 


CHAPTER  H 
ANAEROBIC  DIGESTION 

^  .  1  .   Introduction 

In  order  to  test  the  adaptive  optimization  algorithm,  an  actual 
bioprocesa  was  chosen.  Anaerobic  digestion  was  selected  because  the 
unique  complexities  of  the  process  present  an  interesting  challenge. 

Probably  one  of  the  oldest  uses  of  microbiology,  anaerobic  diges- 
tion converts  a  wide  variety  of  feedstocks  to  a  readily  usable  energy 
source,  methane.  This  ubiquitous  process  exists  in  many  different 
sizes,  from  small  household  composting  units  to  large-scale  industrial 
waste  conversion  systems. 

The  most  frequent  use  is  found  in  municipal  waste  applications 
where  anaerobic  digestion  is  used  to  treat  domestic  sludge  and  biode- 
gradable garbage,  such  as  food  products  and  paper.  A  trip  to  many  farms 
will  find  the  use  of  digesters  to  dispose  of  cow,  chicken  and  pig  man- 
ures as  well  as  excess  crops.  In  fact,  many  researchers  have  studied 
the  feasibiliity  of  using  crops  specifically  to  produce  energy  via 
digestion  [35,36].  Even  the  degradation  of  some  hazardous  wastes  can  be 
accomplished  through  anaerobic  digestion  [37,8].  Overall,  this  process 
is  extremely  important  for  the  removal  of  energy  from  wastes  to  prevent 
imbalancing  the  energy  equilibrium  of  the  environment. 

Anaerobic  digestion  naturally  occurs  in  soil,  lake  and  ocean  beds, 
and  animal  digestive  tracts.  It  can  generally  be  described  as  a  two- 
stage  mechanism  [39].   Complex  organic  substrates  (e.g.,  carbohydrates, 

lipids,  and  proteins)  are  hydrolyzed  and  fermented  in  the  first  stage  by 
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a  group  of  bacteria,  referred  to  as  acidogens,  into  organic  acids, 
carbon  dioxide,  and  hydrogen  gas.  These  products  are  converted  in  the 
second  stage  to  methane  and  carbon  dioxide  by  a  unique  group  of  bacteria 
called  methanogens. 

However,  this  simple  representation  does  not  encompass  all  reac- 
tions occurring  in  a  typical  digester.  For  example,  a  group  of  homo- 
acetogenic  bacteria  reduce  carbon  dioxide  to  acetate  and  accounts  for  an 
appreciable  amount  of  the  methane  produced  [40].  In  addition,  this  two 
stage  description  obscures  the  delicate  interaction  between  acidogenic 
and  methanogenic  bacteria.  Experiments  ran  with  separated  stages  pro- 
duce different  behavior  and  intermediates  [41,42].  This  interaction  is 
also  the  cause  of  many  of  the  problems  associated  with  digester  opera- 
tion. Overfeeding  and  inhibitors,  such  as  oxygen  and  chloroform,  can 
quickly  cause  a  healthy  reactor  to  turn  sour. 

The  health  of  the  anaerobic  digester  is  usually  indicated  in  prac- 
tice by  measurements  of  overall  gas  production,  methane  percentage  of 
effluent  gas,  and  volatile  fatty  acids  (VFAs).  Volatile  fatty  acids  are 
intermediates  and  are  considered  to  be  a  good  indicator  of  imbalance 
[43].  Acetic,  propionic,  n-  and  iso-butyric,  and  n-  and  iso-valeric 
acids  are  usually  measured  with  acetate  usually  being  the  most  predomi- 
nant. A  rise  in  VFA  concentration,  ^specially  the  higher  weight  acids, 
is  usually  a  precursor  to  failing  methane  gas  production.  However,  VFA 
analysis  is  time  consuming  and  usually  only  indicates  a  reactor  in  a 
declining  mode.  In  addition,  the  mechanisms  of  VFA  reactions  are  con- 
tinually being  debated  [44,45].  Most  engineers,  from  a  practical  point 
of  view,  are  concerned  mainly  with  the  methane  gas  production.  This 
value  is  easily  measured  on-line  (see  Appendix  C)  which  is  fortunate 
since  methane  gas  is  the  desired  product. 
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Most  efforts  to  optimizR  the  methane  production  have  followed 
traditional  procedures.  Individual  species  of  bacteria  are  isolated  and 
kinetic  data  are  acquired  through  pure  culture  experiments.  This  ap- 
proach obviously  fails  to  include  interaction  between  species.  In 
addition,  there  are  many  bacteria  that  have  not  been  isolated  ['46]. 

Another  procedure  involves  batch  experiments  [47].  Reactors  are 
primed  with  a  starter  culture  and  a  quantity  of  substrate  under  varying 
conditions.  Methane  production  rates  are  then  measured  throughout  the 
viability  of  the  culture.  This  method  accounts  for  interaction  but  does 
not  provide  enough  predictive  information  on  continuous  operation.  The 
variety  of  microbes  in  a  typical  reactor  is  immense  and  relative  popula- 
tion sizes  are  dependent  on  substrate  concentration.  Therefore,  contin- 
uous operation  can  not  be  properly  described  by  batch  experiments. 

The  most  popular  optimization  method  is  continuous  steady-state 
experiments.  Step  changes  are  made  in  a  particular  manipulated  variable 
and  the  subsequent  steady-state  output  is  measured.  However,  typical 
residence  times  of  10  to  20  days  would  result  in  months  of  waiting  for 
each  steady  state.  During  this  period  of  transition,  many  upsets  could 
occur,  such  as  population  shifts,  to  render  these  results  suspect. 

In  view  of  these  problems  and  complexities,  adaptive  optimization 
could  provide  a  useful  alternative  for  maximizing  production  of  an 
anaerobic  digester.  Dynamic  data  would  be  used  to  predict  steady-state 
behavior,  thus  circumventing  the  need  for  steady  states.  Also,  on-line 
implementation  would  avoid  the  problems  associated  with  ad  hoc  experi- 
mentation. 
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H.2.      Temperature  Optimizablon 

Temperature  is  one  of  the  moat  Important  factors  in  anaerobic 
digestion  that  influences  methane  gas  production.  Enzyme  kinetics, 
dissociation  constants,  and  death  rates  are  greatly  affected  by  small 
changes  in  the  culture  temperature.  Johnson  et  al.  [3^]  proposed  that 
bacteria  are  mostly  influenced  by  changes  in  enzyme  kinetics.  As  tem- 
perature rises,  enzyme  activation  Increases,  while  at  the  same  time 
enzyme  denaturation  increases  (see  Chapter  3).  In  addition,  higher 
temperatures  also  increase  irreversible  destruction  of  these  vital 
proteins.  These  mechanisms  will  cause  a  typical  bacteria  to  have  a 
range  of  viability  as  well  as  an  optimum  temperature  for  growth. 

There  are  two  operating  ranges  that  have  distinct  characteristics, 
thermophilic  and  mesophilic.  Thermophilic  operation  is  usually  between 
50  and  60°C  whereas  mesophilic  digesters  usually  run  between  30  and 
40°C.  Many  authors  [48,49,50]  have  discussed  the  relative  merits  of 
each  range,  concluding  that  thermophilic  digesters  have  the  most  advan- 
tages (e.g.,  better  dewatering  characteristics,  higher  fraction  of 
pathogens  destroyed,  lower  cell  mass  production,  and  increased  rate  of 
digestion)  but  require  stricter  controls  on  temperature  variations. 

A  comprehensive  work  by  Zeikus  [46]  on  the  biological  morphology  of 
methanogens  describes  microbial  species  that  are  specific  to  each 
range.  Many  researchers  [50,51]  have  verified  that  separate  populations 
do  exist  in  actual  digesters.  However,  enumeration  experiments  on  a 
mesophilic  sewage  fermentor  established  that  10  %  of  the  total  bacteria 
population  were  typical  thermophilic  organisms  [50]. 

Much  research  effort  has  been  expended  to  ascertain  the  temperature 
profiles  of  anaerobic  digesters.   Pfeffer  [52]  showed  the  existence  of 
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two  different  optima  for  a  digester  fed  with  shredded  municipal  waste. 
Mesophilic  operation  was  optimized  between  40  and  42°C  and  had  a  higher 
maximum  specific  growth  rate  than  thermophilic  operation,  which  was 
optimized  between  60  and  62°C.  Other  authors,  however,  presented  re- 
sults demonstrating  little  difference  in  gas  production  between  the  two 
ranges  [53,5^,55]. 

Chen  and  Hashimoto  [56]  used  the  data  of  Pfeffer  and  others  [57,58] 
to  propose  the  existence  of  two  ranges.  However,  in  a  later  paper  Chen, 
Hashimoto  and  Varel  [59]  presented  data  on  a  series  of  exhaustive  exper- 
iments to  determine  the  effect  of  dilution  rate  and  temperature  on  beef- 
cattle  manure  (Figure  6-1).  In  this  paper,  they  proposed  that  a  smooth 
transition  actually  existed  between  mesophilic  and  thermophilic  condi- 
tions. Longer  acclimation  periods  were  shown  to  improve  results.  They 
concluded  that  gas  production  increased  steadily  from  low  to  high  tem- 
peratures with  maximization  occurring  between  55  to  60°C. 

These  authors  also  reported  sudden  drops  in  steady-state  gas  pro- 
duction at  particular  temperatures.  More  specifically,  gas  production 
seemed  to  occur  only  for  a  particular  range  with  very  little  change 
within  this  range  (Figure  6-1).  The  range  also  became  more  narrow  as 
the  dilution  rate  increased.  (A  possible  explanation  of  this  phenomenon 
is  proposed  in  Chapter  6.  In  addition,  a  method  of  improving  optimi- 
zation and  avoiding  such  behavior  is  presented.) 

Much  debate  still  presides  over  this  important  question  of  tempera- 
ture profiles.  Adaptive  optimization  offers  an  unbiased  method  of 
providing  the  answer.  In  truth,  much  of  the  data  presented  in  the 
research  described  above  is  specific  to  a  particular  system  and  feed 
substrate.   The  general  approach  of  adaptive  optimization  is  well-suited 
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for  many  systems.   In  Chapter  7,  results  of  this  method  applied  to  a 
continuous  glucose-fed  anaerobic  digester  are  presented. 


CHAPTER  5 
ALTERNATE  SAMPLING  BASIS 

5.1.  General  Description 
Before  a  discrete  model  can  be  implemented  on  line,  a  sampling 
interval  must  be  chosen.  This  value  is  important  to  the  success  of  the 
control  or  optimization  algorithm,  but  is  usually  chosen  haphazardly. 
If  the  sampling  interval  is  too  small,  changes  in  the  measurement  output 
could  be  masked  by  process  noise.  However,  dynamic  information  could  be 
lost  if  too  large  a  value  is  chosen.  Smith  and  Corripio  [16]  have  shown 
that  best  results  can  be  obtained  by  choosing  a  sampling  interval  to  be 
about  ten  to  twenty  percent  of  the  dominant  time  constant. 

A  major  problem  with  this  method  is  the  variability  of  the  time 
constant  as  the  process  changes  from  one  steady  state  to  another.  This 
is  especially  magnified  in  bioprocess  optimization  where  the  bioreactor 
changes  from  a  low  production  state  characterized  by  slow  dynamics  to  a 
high  production  state  with  substantially  faster  dynamics. 

This  problem  may  be  circumvented  by  considering  an  alternate  sam- 
pling basis,  i.e.,  transforming  the  independent  variable  from  time  to 
some  other  variable,  <J).   This  new  basis  can  reduce  the  variation  of  the 

time  constant. 

This  transformation  can  be  viewed  in  general  terms  by  considering 
the  following  general  single-input  single-output  system: 


dx  „,  s 
-rr  =  f(x,u) 
dt 
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where  x  is  the  output  measurement  and  u  is  the  control  input.   Using  a 
first  order  Taylor  series,  the  time  constant,  t,  can  be  obtained  from: 


3f 
s 

3x 


]  = 


-1 


Here,  \   is  the  eigenvalue.   The  transformation  is  performed  by  dividing 
f  by  the  derivative  of  the  new  basis  with  respect  to  time. 


dx 


r  d(J) ' 


•1 


-  =  g(x.u)  =  f[^j 


The  new  eigenvalue  can  be  easily  related  to  the  time  based  eigenvalue. 


dt 


df 

■ 3 

•  3x 


Assuming  the  transformation  is  performed  for  a  stable  system,  the  fol- 
lowing condition  must  be  met  in  order  to  maintain  stability. 


This  condition  implies  that  the  new  eigenvalue  is  always  negative. 
Therefore,  the  new  independent  variable  must  be  a  monotonically  increas- 
ing function  of  time.  This  can  be  more  clearly  understood  by  consider- 
ing the  following  example. 


5.2.   Example 
A  constant  volume  isothermal  continuous  stirred  tank  reactor  (CSTR) 
in  which  an  n    order  reaction  takes  place  can  be  described  by  the 
following  continuous-time  dynamic  model: 
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V  ^  =  FC   -  FC  -  Vkc"  (5-1  ) 

dt     o 


where  C  is  the  reactant  concentration,  C^  the  feed  reactant  concentra- 
tion, V  the  reactor  volume,  F  the  flow  rate  and  k  the  kinetic  con- 
stant. Let  C  be  the  measured  variable  and  F  the  manipulated  variable, 
We  can  change  equation  (5-1)  to  dimensionless  form  by  defining 


c  =  ^  (5-2) 

C 
o 


f  =  ^  (5-3) 

o 


6  =  tkC^"^  (5-4) 

o 


in  which  case  equation  (5-1)  becomes; 


^  =  f  -  fc  -  c"  (5-5) 

d9 


At  steady  state  the  above  equation  yields 


n 
c 

f  =    " 


S     1  -  G 

s 


(5-6) 


where  the  subscript  s  denotes  steady  state.   The  eigenvalue  is 


A  =  -(f   +  nc""^)  (5-7) 

s     s 


which  corresponds  to  the  time  constant 
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1  ^  "  '^ 

T.  =  — —T   = ^— T      (5-8) 

1    -  ^   n-1  ,,   ,  n     n-1 

f  +  nc  (1 -n)p  +  nc 

S3  3      3 

This  is  clearly  a  function  of  the  concentration  at  the  operating  steady 
state  c„  =  (0,1).  For  example,  if  n  =  3,  for  ^%  conversion  (c^  =  .99) 
the  time  constant  is  0.01  while  for  99?  conversion  (c  =  .01)  it  in- 
creases to  3320.  The  lower  the  concentration,  the  lower  is  the  reaction 
rate  and  the  slower  are  the  dynamics  as  evidenced  by  the  substantial 
increase  in  the  time  constant. 

Let  6,(6)  be  the  ( dimensionless)  cumialative  amount  of  reactant 
consumed  until  time  9.   Then 


J  (6)  =  /  c"  de  (5-9) 

0 


and 


^=c"  (5-10) 


Clearly  B  is  a  monotonically  increasing  function  of  6.   3y  dividing 
equation  (5-5)  by  equation  (5-10)  we  obtain 


dc  ^  fc""  -fc'^"-^^  -  1         (5-11) 
dS^ 


What  we  have  done  is  change  the  independent  variable  from  time  to  cumu- 
lative amount  of  reactant  consumed.   The  eigenvalue  of  (5-11)  is 


f 

X^  =  -  -^   [n  ^  (l-n)c^]        (5-12) 

c 
s 
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and  using  equation  (5-6)  we  obtain  the  following  "time"  constant; 


1    G  (1-0^) 

1  s s 

"^2  "  ~  IT  "  n  +  (1-n)c 

2  3 


(5-13) 


This  also  varies  with  c^,  but  considerably  less  than  x,.  For  example, 
if  n=3,  in  the  range  c^  =  [0.01,  0.99]  the  maximum  value  is  0.134  and 
the  minimum  is  0.0033,  a  ratio  of  '40   versus  a  ratio  of  332,000  for  t^. 

Let  6-,(9)  be  the  (dimensionless)  cumulative  net  inpuc  (input  - 
output)  of  reactant,  i.e. 


(9)  =  J   f(1-c) 
*       0 


de 


(5-14) 


This  is  also  a  monotonically  increasing  function  of  e.   Dividing  equa- 

dB 
tion  (5-5)  by  —^   we  obtain: 


do 
d6. 


f(1-c) 


(5-15) 


The  "time"  constant  is 


c  (1-c  ) 
s    s 

n  +  (1 -n)c 


(5-16) 


i.e.,  it  is  equal  to  j    , 

Finally,      let   S;,(9)    be     the     (dimensionless)     cumulative     amount     of 
reactant   fed,    i.e. 


(e)  =  /    f-1  de 

0 


(5-17) 
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Again  we  have  a  monotonioally  increasing  function  of  G,  and  when  we 
change  from  time  to  6^  as  the  independent  variable,  equation  (5-5) 


becomes 


^'^  =  1  -  c  -  ^  (5-18) 


dBij  f 


and  the  corresponding  "time"  constant  is 


G 

s 


^4  =  n.  (1-n)n  ^5-^9) 


Table  5-1  presents  the  ratio  of  the  maximum  "time"  constant  to  the 
minimum  "time"  constant  versus  n  for  c,  =  [0.01,  0.99].  It  makes  it 
clear  that  by  switching  the  independent  model  variable  from  time  to 
cumulative  amount  of  reactant  consumed,  cumulative  net  amount  of  reac- 
tant  fed  or  even  cumulative  amount  of  reactant  fed,  the  variability  of 
the  "time  constant"  with  c  is  considerably  reduced.  This  is  increas- 
ingly the  case  as  the  reaction  order,  and  with  it  the  nonlinearity  of 
the  process,  increases. 


Table  5-1 .   Ratio  of  maximum  to  minimum  "time"  constants 
versus  reaction  order 


n 

T. ,  max 
T  ,  min 

'2' 
'2' 

max 
rain 

T^,  max 
T-,  rain 

T  j^ ,  ma  X 
T^^,  min 

1 

99 

25.3 

99 

2 

4974 

34.5 

195 

3 

332088 

40.3 

289 

H 

2.5  X  lo"^ 

44.6 

382 
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The  above  observation  motivates  discretization  on  the  basis  of 
model  (5-11)  instead  of  model  (5-5).  Using  the  Euler  method  an  approxi- 
mate discretization  of  (5-11)  is: 


ciQ^+^&^)   =   0(3^)  + 


A62  [f(62)c(62)""  -  f (22)0(32)"^'''"^^  -  1]  (5-20) 


If  a  control  or  optimization  algorithm  employs  the  above  model,  sampling 
occurs  when  a  fixed  amount  of  &  has  accumulated  (constant  A32)  instead 
of  when  a  fixed  amount  of  time  has  passed  (constant  A9).  This  in  effect 
creates  a  variable  sampling  period.  The  sampling  interval  is  automati- 
cally lengthened  or  shortened  as  dictated  by  the  change  in  the  speed  of 
the  process  dynamics. 

This  concept  is  of  course  not  restricted  to  the  above  mentioned 
reactor  example.  For  many  processes  which  exhibit  wide  variations  of 
the  effective  time  constant,  the  associated  problems  can  be  reduced  by 
considering  as  the  independent  variable  the  cumulative  amount  of  a 
quantity  generated,  consumed  or  fed,  rather  than  time.  In  this  manner 
the  sampling  period  will  be  automatically  adjusted,  increasing  when  the 
dynamics  slow  down  and  decreasing  when  they  speed  up. 

5.3.  Application  to  Anaerobic  Digestion 
The  anaerobic  digestion  process,  in  which  cellulosic  material  or 
sugars  are  converted  to  methane  and  carbon  dioxide,  can  exhibit  consi- 
derable variation  in  the  effective  time  constant.  The  lower  the  dilu- 
tion rate  (flow  rate /volume)  is,  the  lower  is  the  methane  production 
rate  and  the  slower  are  the  dynamics. 
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To  exemplify  the  advantages  of  alternate  sampling  basis  for  anaero- 
bic digestion,  an  adaptive  control  simulation  is  presented  in  this 
section.  Many  excellent  studies  have  been  published  for  adaptive  con- 
trol of  biochemical  systems  [60, 6l]  as  well  as  anaerobic  digestion 
[62,63].  The  adaptive  control  law  used  here  is  a  simple  quadratic  type 
since  the  main  purpose  is  to  exhibit  the  properties  of  an  alternate 
sampling  basis. 

Whether  the  purpose  of  the  digester  operation  is  energy  production 
or  wastewater  treatment  or  both,  a  reasonable  control  objective  is  to 
attempt  to  maintain  a  specified  methane  production  rate.  (A  word  of 
caution  though:  If  the  digester  is  subject  to  receiving  inhibitory 
material  in  the  feed,  maintenance  of  a  constant  methane  production  rate 
should  not  be  the  primary  control  objective  [64]).  The  methane  produc- 
tion rate  can  be  easily  measured  on  line.  A  simple  scheme  is  described 
in  Appendix  C.  An  easy  to  manipulate  variable  that  substantially  af- 
fects the  methane  production  rate  is  the  dilution  rate.  Thus  the  objec- 
tive is  to  control  the  methane  production  rate  by  manipulating  the 
dilution  rate. 

Anaerobic  digesters  have  a  complex  mixed  culture,  and  changes  in 
environmental  conditions  may  cause  population  changes.  Furthermore,  the 
feed  is  often  ill-characterized  and  changing.  For  these  reasons  an 
adaptive  controller  is  employed  instead  of  basing  the  controller  design 
on  a  fixed  process  model.  A  fixed  process  model  [65]  is  used,  though, 
for  simulating  the  process. 

The  adaptive  controller  is  based  on  the  following  on-line  identi- 
fied model: 
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Q^^^(B^-A8)  =  e^QpH^^^''  *  ®2'^^^^  *  ®3  ^^"^^^ 


where  Q^uj^  is  the  methane  production  rate  (liters  of  methane  per  day),  D 
is  the  dilution  rate  and  6  are  the  on-line  estimated  parameters. 
Instead  of  time,  the  independent  variable  g  is  the  cumulative  amount  of 
methane  produced,  something  easily  measured  on  line.  Sampling  takes 
place  whenever  B  increases  by  A6. 

The  adaptive  controller  employed  is  a  quadratic  self-tuning  con- 
troller [66]  based  on  the  minimization  of  the  performance  measure: 


^^°^»^^  =  ^QCH4,  desired  "  ^CH^^^*^^^^' 


+  q[D(B)  -  D(B-A6)]^  (5-22) 


The  resulting  control  law  is: 


D(B)  =  ^  {qD(B-AB)  + 

q  +  62 


where  6.  denotes  the  current  parameter  estimate. 

The  parameter  estimator  is  a  modification  of  the  recursive  least 
squares  with  variable  forgetting  factor  of  Fortescue  et  al.  [22].  The 
modification  is  that  the  m  (m  ^  number  of  parameters)  most  recent  data 
are  excluded  from  forgetting.  This  has  the  advantage  that  even  at 
periods  of  high  forgetting,  the  estimator  retains  the  minimum  amount  of 
information  required  (or,    equivalently,  the  estimator  can  employ  higher 
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rates  of   forgetting  of  past    data).      A  recursive   formulation   that  accom- 
plishes   this    is   the    following: 

At  each  sampling   instant  v: 


T  * 

e         =  y         -  X  9 

v-m         v-m       — v-m  - 


P       X 

K  • 


^-"^        1    -    x^        P*x 

—v-m       -v-m 


»  * 

6      '-   6      +   K  e 

—         —  v-m     v-m 


*  T  * 

P     -^    [I   -  K  X        ]P 

v-m  —v-m 


e        =9 

-v-m       - 


* 

P  =   P 

v-m 


For   i   =    1    to   m 


T 
e  .    =   y(v-m+i)   -  x  .9        ^,    . 

v-m+i  -v-m+:-v-m+i-1 


P  X 

v-m+i -1   —v-m+i 
K 


v-ra+i        1    +   x'^  P    '  x 

-v-m+i     v-m+i-1   -v-m+i 


i  =      ft  +      IC  P 

-v-m  +  i       -v-m+i-1  v-m  +  i     v-m+1 


T 

p         =ri  —  K         X        ]p 

v-m+i  v-m+i   —v-m+i        v-m+1-1 
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Next  i 

2 
e 

»  =  , :^ 


P   -  P  /A 


The  notation  is  as  follows: 


V         =   independent  variable  enumerating  the  number  of  sampling 
periods  (At  for  conventional  algorithms  or,    in  our  case, 
A8) 

_9         =   parameter  estimate  vector  obtained  with  forgetting  and 

without  considering  the  last  m  data 

« 

P         =   corresponding  covariance  matrix 

Yy.jj       =   measurement  at  sampling  instant  v-m 

QQjjji|(3-tnA6)  in  our  case 

X         =   information  vector 
—v-m 

[QQ^i,( v-m-1  ) ,  D(v-m-l),  1]   in  our  case 
e^_jjj       =   estimation  error 

K-m  =       gain 


v-m 


-v-m+1 


update  of  9  using  y(v-m+1)  and  without  forgetting 


9    .     =   update  of  9    .  ,  using  y(v-m+i-l)   and  without  for- 
-v-m+i         '^        — v-ra+i-1     o     j^ 

getting 

^v-m+i     ^   corresponding  covariance  matrix 

9^         =   final  parameter  estimate  at  instant  v  (used  in  the  con- 
trol law) 
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X  =       forgetting  factor.   The  constant  a   adjusts  the  speed  of 

forgetting 

Essentially  what  this  algorithm  does  is  update  the  value  Q     that  the 

parameter  estimates  would  have  if  the  last  m  data  were  not  received, 

using  a  variable  forgetting  factor.  Then,  the  last  m  data  are  used,  one 

at  a  time,  to  update  the  parameter  vector  without  forgetting.   It  should 

be  remarked  that  only  P  ,  Q_     and  the  inputs  and  outputs  contained  in 

X  X    ,  are  stored  after  each  iteration;  so  the  algorithm  requires 

-v'    -v-m+1 

approximately  as  much  storage  as  a  conventional  least  squares  estimator. 
Figures  5-1  to  5-3  show  simulation  results  using  the  Smith  et  al. 
[65]  model  for  a  5  liter  glucose-fed  digester  with  glucose  feed  concen- 
tration of  30  g/L.  The  sampling  period  was  selected  to  be  1  liter  of 
methane  produced.  The  algorithm  parameters  were  m=3,  o=0.005  and 
q=10,000.  The  initial  set-point  was  0.93  L/day  and  after  an  initializa- 
tion period  of  10  sampling  instants  (during  which  the  system  was  sub- 
jected to  random  input  excitation),  it  was  changed  to  6  L/day.  Subse- 
quently, the  set-point  was  returned  to  its  original  value.  As  Figure  5- 
1  shows,  the  controller  successfully  drives  the  methane  production  rate 
to  the  desired  values.  Figure  5-2  shows  the  corresponding  hydraulic 
retention  time  (inverse  of  dilution  rate).  Figure  5-3  exhibits  the 
time-variation  in  the  sampling  period.  Keeping  A6  constant  at  1  L 
results  in  equivalent  time  sampling  period,  At,  varying  from  4  hours 
(when  the  system  dynamics  are  faster  at  the  high  methane  production  rate 
set-point)  to  26  hours  (when  the  dynamics  are  slower  at  the  low  methane 
production  set-point).  Finally,  Figure  5-4  shows  that  if  a  constant  At 
self-tuning  controller  were  used,  with  At  equal  to  that  of  the  constant 
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Figure     5-2.    Hydraulic     retention 
sampling  basis. 
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Figure  5-3.  Time  variation  of    sampling  interval  of  simulation  using 
alternate  sampling  basis. 
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Figure  5-'^.   Gas  production  of  simulation  using  time  sampling  basis, 
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AB  algorithm  in  the  initial  steady  state  (and  identical  algorithm  tun- 
ing), the  system  does  not  speed-up  the  sampling  rate  and  as  a  result 
reaches  set-points  considerably  more  slowly. 

The  advantages  gained  by  this  transformation  are  basically  three- 
fold. Variation  in  parameter  estimates  will  be  reduced  because  of  the 
reduction  in  the  variation  of  the  time  constant.  The  system  will  be 
better  represented  by  a  linear  discrete  model  which  will  lead  to  im- 
proved stability.  Finally,  the  time  required  to  optimize  a  system  will 
be  reduced  because  the  sampling  period  will  get  smaller  as  the  dynamics 
of  the  system  increase. 

It  is  also  necessary,  however,  to  consider  the  potential  problems 
associated  with  the  transformation.  A  reliable  measurement  of  the 
independent  variable  is  required.  For  the  anaerobic  digestion  process, 
this  is  not  a  problem  since  methane  gas  can  be  measured  in  increments 
very  easily  (see  Appendix  C).  Also,  all  changes  to  process  inputs  must 
be  made  with  respect  to  the  new  basis.  For  continuous  anaerobic  digest- 
ers, this  approach  leads  to  an  attractive  operation  mode,  constant 
reactor  yield,  which  is  described  in  the  next  chapter.  Lastly,  because 
the  sampling  interval  is  variable  with  respect  to  time,  it  must  always 
be  large  enough  to  allow  for  data  retrieval  and  manipulation. 


CHAPTER  6 
CONSTANT  YIELD  OPERATION 

The  majority  of  anaerobic  digesters  are  operated  in  a  continuous 

semibatch  mode.   Substrate  solutions  are  usually  prepared  and  fed  on  a 

daily  basis.   The  amount  fed  is  maintained  at  a  constant  level  thus 

keeping  the  hydraulic  residence  time  (HRT)  constant.    For  a  simple 

Ghemostat  type  fermenter  (as  in  Figure  C-1),  the  HRT  is  defined  by  the 

following  relationship: 


HRT  =  V/F  =  D~^  (6-1 ) 


where  V  is  the  liquid  volume  of  the  reactor,  F  is  the  volume  of  feed  per 
time,  and  D  is  the  dilution  rate.  The  kinetics  of  a  continuous  fer- 
menter can  be  described  as: 


^   =  m(s)x  -  Dx  -  k^X  (6-2) 

at  d 


^  ._   -  )^   .  D(s  -S)  (6-3) 

dt      Y         o 


Q   =  Y  u(S)X  (6-4) 

P    P 


Here,  X  is  the  microbial  concentration,  p  is  the  specific  growth  rate, 

kj  is  the  death  rate,  S  is  the  substrate  concentration,  Sq  ^^  ^^^   inlet 

concentration  of  substrate,  Y^^  is  the  biomass  yield  (mass  of  cells 

formed/mass  of  substrate  consumed),  Qp  is  the  product  gas  rate,  and  Yp 
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is  the  product  yield  (volume  of  gas/mass  of  cells  formed).  In  many- 
digestion  systems,  methanogenesis  can  be  assumed  to  be  the  rate  limiting 
step.  Therefore,  X  is  the  methanogenio  bacteria  concentration  (g/L),  S 
is  the  total  volatile  fatty  acid  concentration  (g/L  COD),  and  P  is  the 
methane  (L). 

As  discussed  in  Chapter  H,  digester  failure  is  mostly  caused  by 
overfeeding.  An  upset  in  the  reactor  lowering  the  specific  growth  rate 
will  result  in  rising  VFA  concentrations  and  a  declining  methanogenio 
bacteria  concentration.  Furthermore,  substrate  inhibition  may  cause  the 
specific  growth  rate  to  decrease  even  more.  Consequently,  a  digester 
may  appear  healthy  but  in  actuality  be  close  to  souring.  To  prevent 
failure,  digesters  are  operated  at  sub-optimal  HRTs  and  methane  produc- 
tion is  watched  very  closely. 

This  behavior  has  also  been  well  documented  in  research  lab  digest- 
ers. Investigators  have  noticed  slight  changes  in  a  controlled  vari- 
able, such  as  temperature  or  pH,  will  cause  a  sudden  collapse  in  gas 
production  [^3,53].  For  example,  in  the  study  provided  by  Chen  et  al. 
[59],  digesters  were  operated  at  different  temperatures  (see  Figure  6- 
1).  These  researchers  noticed  that  steady-state  gas  production  was 
minimally  affected  within  a  certain  range  of  temperatures  for  a  given 
HRT.  Outside  this  range  gas  production  was  practically  zero.  They  also 
noted  that  this  range  narrowed  as  the  HRT  decreased. 

These  phenomena  can  be  predicted  by  substituting  a  simple  substrate 
inhibition  functionality  for  the  specific  growth  rate: 


y(S)  =  ^H! (6-5) 

1    s   S 
1  +  —  +  — 

S   K 
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Figure    6-1  . 


Experimental  results  of  [59].  Temperature  effects  on  the 
steady-state  methane  production  for  different  retention 
times. 
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where  K  and  Kx  are  kinetic  parameters,  and  using  the  thermodynamic 
functionality  of  equation  (3-24).  Temperature  profiles  of  steady-state 
gas  production  predicted  from  this  model  are  shown  in  Figure  (6-2). 
Although  these  process  curves  violate  common  logic,  an  explanation  can 
be  given  by  considering  the  effect  of  temperature  and  substrate  concen- 
tration on  the  specific  growth  rate.   If  the  temperature  is  increased 

toward  the  optimum,  p  of  equation  (3-24)  will  increase.   However,  the 

m 

substrate  level  will  decrease  so  as  to  maintain  the  product  term  of 
equation  (6-4)  at  an  approximately  constant  value.  At  the  drop  off 
points  away  from  the  optimum,  rising  substrate  levels  due  to  a  decrease 

in  M 

m 

begin  to  have  a  dramatic  inhibitory  effect  on  the  growth  rate.  Subse- 
quently, faster  feeding  rates  (lower  HRTs)  will  amplify  this  effect  and 
result  in  a  reduced  range  of  possible  operating  temperatures.  One  can 
easily  imagine  the  problems  of  trying  to  operate  a  digester  at  the 
maximum  feed  rate  (lowest  HRT). 

A  way  of  operating  the  digester  that  could  possibly  avoid  these 
problems  is  the  constant  reactor  yield  feeding  mode.  In  this  method, 
the  amount  of  substrate  fed  per  amount  of  methane  produced  is  held 
constant  (see  Appendix  C  for  implementation).  This  functionality  is 
given  by  the  following  equation: 

Q     Y  p(S)X 
-2-  =  _P (6-6) 


DS  V    DS  V 
o       o 


where  Y  is  the  reactor  yield.  Substituting  into  equations  (6-2)  to  (6-4) 
gives 
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Figure  6-2.  Theoretical  effect  of  temperature  on  steady-state  gas 
production  of  a  constant  HRT  digester  (at  different  HRTs) 
using  equations  (3-25),  (3-26),  (6-2),  (6-3),  (6-4)  as 
process  model. 
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f  =  M(S)X  (1  -:^X)  -K^X  (6-7) 

o 

(S  -S)Y  Y 
dS  ^  ia(S)X  r      o   '       X  p  _  1  ,._g. 

dt    Y     '■   S   YV       -' 
X       o 


The  steady-state  values  of  X,  S  and  methane  rate  are; 


P 

YVS 
X  p 


«PS=  ^^^^s^  -  ^d^^V  ^^-"^ 


Equation  (6-10)  can  also  be  rearranged  in  convenient  terms  of  substrate 
conversion  fraction. 


—2 =  -TL  (6-12) 

S      Y  Y  . 

X  p 


Constant  yield  operation  has  the  advantage  that  the  substrate  concentra- 
tion is  held  constant  thus  avoiding  substrate  inhibition  effects.  In 
addition,  for  a  small  death  rate,  the  methanogenic  bacteria  concentra- 
tion is  also  constant  and  therefore  Q^  is  affected  by  y  only.   This 

p  m 

results  in  a  smooth  parabolic  process  curve  with  a  clearly  distinguish- 
able optimum. 

For  the  above  case  of  temperature  optimization,  this  characteristic 
has  a  dramatic  effect.  In  Figure  (3-3)  digester  temperature  profiles 
are  shown  using  constant  yield.  The  advantage  gained  here  is  that  the 
optimum  is  more  pronounced  and  easier  to  locate.   Also,  the  digester  can 
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Figure  5-3.  Theoretical  effect  of  temperature  on  steady-state  gas 
production  of  a  constant  yield  digester  (at  different 
yields)  using  equations  (3-25),  (3-26),  (o-2),  (6-3),  (6-4) 
as   process  model. 
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be  operated  at  high  production  rates  with  less  chance  of  sudden  fail- 
ures. 

The  advantages  of  constant  yield  operation  are  especially  important 
for  the  case  of  inhibitors  in  the  feed.  Harmon  et  al.  [67]  have  shown 
results  on  a  simulated  anaerobic  digestion  system  [55]  where  an  inhibi- 
tor was  input  into  a  reactor  through  the  feed.  The  constant  HRT  reactor 
ultimately  failed  whereas  the  constant  yield  reactor  had  only  a  slight 
rise  in  VFA  concentration.  In  this  practical  case,  constant  yield 
provides  an  automatic  method  of  raising  the  HRT  for  declining  gas  pro- 
duction. 

It  should  be  noted  that  the  above  results  assume  the  inlet  sub- 
strate concentration  remains  constant.  In  actual  systems,  feedstocks 
may  vary  and  thus  produce  negative  results  under  constant  yield  opera- 
tion. For  example,  an  increase  in  substrate  concentration  will  result 
in  a  rise  in  gas  production  which  will  also  increase  the  feed  rate. 
Consequently,  constant  yield  policy  will  overfeed  the  reactor  [64]. 
However,  Pullammanappalli 1  et  al.  [64]  have  presented  an  expert  system 
which  can  handle  such  upsets. 


CHAPTER  7 
EXPERIMENTAL  RESULTS 

7.1.   Inbroductlon 

In  the  experiments  performed  on  the  anaerobic  digester  described  in 
Appendix  C,  steady  state  methane  production  rate  was  the  performance 
index  and  temperature  was  the  control  variable.  The  digester  was  oper- 
ated on  a  constant  yield  basis  (see  Chapter  6)  and  allowed  to  come  to 
steady  state  before  the  start  of  the  experiment.  VFAs  and  methane  gas 
rate  were  constant  for  ^48  hours  to  verify  the  steady  state. 

Both  the  feed  pump  and  gas  meter  were  calibrated  prior  to  the 
start.   The  yield  is  calculated  from  these  values. 


y  ^   g  (mis  of  methane)  ,   , 

Nr^   (mis  of  feed)  . 


where  r  is  the  meter  calibration  (mis  methane  per  meter  reset),  r^  is 
the  pump  calibration  (mis  feed  per  revolution),  and  N  is  the  number  of 
revolutions  of  the  feed  pump  per  meter  reset.  The  constant  yield  feed- 
ing policy  is  implemented  by  turning  the  feed  pump  on  for  N  revolutions 
each  time  the  gas  meter  resets. 

The  reactor  yield  can  be  transformed  into  chemical  oxygen  demand 
(COD)  conversion  percentage  by  the  following  equation: 


CODJ  =  filter  Of  feed^.    g  COD      ^^_^^^  ^,_,^ 

38.9  g  COD    0.35  liter  CHj^-' 
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Here,  the  first  bracketed  term  is  the  COD  feed  concentration  determined 
in  Appendix  C,  the  second  term  is  the  theoretical  COD  conversion  of 
methane  at  273  K  and  1  atmosphere,  and  the  last  term,  where  T  is  the 
meter  temperature,  corrects  the  gas  volume  measurement  for  temperature. 

The  algorithm  described  in  section  7.2  was  implemented  using  the 
alternate  sampling  basis  described  in  Chapter  5.  The  sampling  inter- 
val, 5,  is  given  by 


<5  =  Zr   (mis  of  methane)  (7-3) 


where  Z  is  the  number  of  meter  resets  per  sampling  period.  At  the  end 
of  each  period,  the  instantaneous  methane  rate,  m,  was  calculated  by  a 
weighted  average  of  the  last  ten  gas  rate  measurements 


L      X^~J  m. 
i  =  1 


where  m.  is  the  methane  rate  of  the  ith  meter  reset,  At.  is  the  time 
between  the  j  and  j-1  resets,  and  w  is  the  weight  (w=0.9  throughout  all 
experiments).   A  weighted  variance  was  also  calculated  to  remove  errors. 


P    10   .        z 

S-"  =  (  Z  A-^-l]"^  [  z   A^"-^  (ra.-m)'^]        (7-5) 
j=1         z-9       ^ 


If  any  measurements  were  more  than  one  standard  deviation  from  m,  they 
were  rejected  and  equation  (7-^^)  was  recalculated. 
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The  temperature,  T,  for  the  sampling  period  was  calculated  by  a 
simple  average  of  the  temperatures  measured  at  each  meter  reset  (see 
Appendix  C)  after  skipping  the  first  ten. 


T  =      I        T.  (7-6) 

i  =  11   ^ 


The  hydraulic  residence  time  was  calculated  by 

z 

V   E   At. 


«^^  =  -zri^-  (7-7) 


where  V  is  the  liquid  volume  of  the  reactor. 

7.2.   Algorithm 
The  model  equation  for  the  experimental  runs  was 


y(i)  =  ay(i-l)  +  bu(i-1 )  +  cu~(i-l  )  +  d  =  Q^h,  _^  (7-8) 


where  y(i)  is  the  gas  measurement  from  equation  (7-*^)  scaled  with  the 
starting  steady-state  gas  rate  and  u(i)  is  the  temperature  average  from 
equation  (7-6)  scaled  with  the  starting  temperature.  At  each  sample, 
the  algorithm  receives  a  new  measurement  and  estimates  new  parameters 
using  the  recursive  least  squares  routine  detailed  in  Chapter  5.  This 
estimator  is  similar  to  the  method  proposed  by  Fortesque  et  al.  [22]  but 
includes  two  modifications.  The  first  involves  resetting  the  covariance 
matrix  when  the  a  posteriori  error  becomes  significantly  high  (see 
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section  3.2).  Resetting  allows  for  faster  discarding  of  data  however 
when  it  is  applied  to  the  most  recent  covariance  matrix  it  discards  all 
data.  This  results  in  periods  where  parameter  estimates  will  be  poor 
until  sufficient  data  is  retrieved.  To  avoid  this  problem,  the  second 
modification  of  applying  the  forgetting  factor  to  a  previous  covariance 
matrix  is  included  (see  Chapter  5).  In  this  alteration,  the  n  most 
recent  data,  where  n  is  the  number  of  parameters  to  be  estimated,  are 
not  discarded. 

The  step  size  variable,  h,  from  Chapter  3  is  used  to  determine  the 
confidence  in  the  new  estimated  parameters.  This  variable  is  calculated 
from  equations  (3-17)  and  (3-18)  with  two  additions.  The  saturation 
parameter  of  equation  (3-17),  ct,  was  constant  as  opposed  to  the  variable 
functionality  of  Appendix  A.  This  change  was  imposed  to  simplify  the 
algorithm.  In  addition,  the  filtered  a  posteriori  error,  v,  was  not 
allowed  to  be  greater  than  100  a.  Errors  larger  than  this  value  would 
cause  h  to  be  1  for  extended  periods. 

As  the  a  posteriori  error  drops,  h  will  approach  0  indicating 
increasing  confidence  in  the  predicted  model.  If  h  drops  below  HOFF, 
the  model  parameters  are  considered  to  have  converged  and  updating  is 
stopped  (i.e.,  the  parameter  estimator  is  turned  off  until  h  rises  above 
HOFF). 

The  predicted  optimum  temperature  is  calculated  using  the  new 
parameters  from  equation  (3-1^). 


T  ,  =  -  b/2c  (7-9) 

opt 
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Here,    "    indicates    estimated    values.      Likewise,    the   necessary    condition 
for  optiraality  from  equation    (3-23)    is 


-^  <    0  (7-10) 


1-a 


After  calculating  the  above  values,  the  algorithm  selects  a  new  set 
point  temperature.  In  simulations  presented  in  Chapter  3,  this  was  a 
relatively  simple  task  but  an  actual  real-time  process  requires  a  much 
more  complicated  scheme.  Certain  limits  must  be  placed  on  the  changing 
of  the  temperature  set  point. 

A  minimum  movement  limit  was  placed  on  the  temperature  set  point 
change  for  two  reasons.  First,  the  temperature  change  cannot  be  less 
than  the  accuracy  of  the  controller.  Second,  the  set  point  change 
should  be  large  enough  to  produce  a  response  in  the  gas  measurement  that 
will  not  be  masked  by  process  noise.  This  limit  was  determined  to  be 
0.5°C.  The  resolution  of  the  controller  further  constrained  the  preci- 
sion to  0.1 °C.  In  other  words,  temperature  had  to  be  set  in  increments 
of  0.1°C. 

The  sensitivity  of  the  process  required  a  limit  to  be  placed  on  the 
maximum  allowable  set  point  change.  Large  changes  in  temperature  could 
shock  the  culture  and  produce  an  unfavorable  response  [49].  This  value 
was  set  at  2°C. 

In  addition,  these  bounds  made  the  implementing  of  the  (^ontrol  law 
of  equation  (3-16)  impractical.  A  new  strategy  was  developed  using  four 
control  modes  to  adjust  the  temperature  set  point.  A  schematic  is  shown 
in  Figure  7-1 . 
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path1-lfh<H2  and  lal  <1 

path2-lfh>H2  or  lal    >  1 

path  3  -  If  h  <  HOFF  and   T  „    =  T   , 

base        opt 

path  4  -  If  h  >  HOFF 


Figure     7-1.    Schematic     diagram     of    e 
algorithm. 


xperimental    temperature    set    point 
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The  first  mode,  INIT,  involves  initialization  and  is  used  only  in 
the  first  six  Intervals  of  the  optimization  run.  In  this  section  of  the 
algorithm,  the  temperature  is  randomly  excited  about  a  base  temperature. 


T   (i)  =  T^    +  NSEQ(i)  (7-11 ) 

sp       base 


Here,  T  is  the  setpoint  at  sample  i,  T^^^gg  is  the  base  (starting) 
temperature,  and  NSEQ  is  the  ith  member  of  a  discrete  white  noise  se- 
quence (zero  mean,  0.5°C  standard  deviation).  Basically,  this  section 
is  used  to  acquire  reasonable  estimates  before  attempting  to  locate  the 
optimum. 

After  initialization  the  algorithm  moves  form  INIT  to  the  identifi- 
cation section,  IDENT.  The  control  law  of  this  mode  is  the  same  as 
equation  (7-11)  and  is  used  when  confidence  in  parameter  estimation  is 
low.  The  step  size  variable,  h,  described  above  is  utilized  as  an 
indication  of  estimation  confidence.  In  other  words  when  h  is  close  to 
1 ,  temperature  set  point  is  adjusted  in  the  IDENT  mode. 

When  h  drops  below  a  predetermined  value,  HI ,  and  the  model  equa- 
tion is  stable  (i.e.,  |a|  <  1),  the  algorithm  uses  the  MOVE  mode  to 
change  the  reactor  temperature.  This  is  accomplished  by  changing  T^^^g^ 
in  equation  (7-11)  and  reducing  the  maximum  noise  to  +/-  0.5°C.  The 
magnitude  of  this  change  further  depends  upon  h.  If  h  is  less  than  H2 
then  Tj33gg  is  changed  by  2°C,  otherwise  the  change  is  1°C.  The  direc- 
tion of  change  is  determined  by  the  optimality  condition  (7-10),  If  the 
second  derivative  is  less  than  0,  T^^^g^  is  moved  toward  the  calculated 
optimum  of  equation  (7-9).  For  a  positive  second  derivative,  equation 
(7-9)  is  the  predicted  minimum  and  the  move  is  made  away  from  this 
value. 
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From  the  MOVE  mode,  the  algorithm  can  transfer  control  back  to  the 
IDENT  mode  in  the  event  that  h  >  1  or  a  >  1  .  However,  if  h  drops  below 
the  estimator  cutoff  value,  HOFF,  and  Tj^^g^  has  converged  to  the  pre- 
dicted optimum,  the  algorithm  removes  the  random  excitation  and  main- 
tains the  temperature  at  the  predicted  optimum.  This  mode  is  called 
CONV  and  is  used  to  maintain  the  temperature  and  monitor  h.  Should  h 
rise  above  HOFF,  set  point  adjustment  will  switch  to  the  IDENT  mode. 

The  algorithm  can  be  summarized  in  the  following  steps: 

STEP  1  :  Obtain  new  gas  measurement  and  temperature  set  point  using 
equations  iJ-H)    to  (7-6) 

STEP  2:   Use  Estimator  of  Chapter  5  to  estimate  new  parameters 

STEP  3:  Calculate  h  form  equations  (3-17)  to  (3-19).  If  h  <  HOFF  for 
three  sampling  intervals  do  not  update  parameters  and  for- 
getting factor  is  set  to  one  for  next  sample 

STEP  4:   Calculate  new  set  point  using  equation  (7-11)  and  Figure  7-1 

STEP  5:   Go  to  STEP  1 

Algorithm  parameters  are   given   in   Table   7-1. 

7.3.  Tnermopnilic 
The  first  run  was  performed  at  a  starting  temperature  of  49°C.  The 
reactor  yield  was  11.03  liters  of  methane  per  liter  of  feed.  COD  con- 
version percentage  (determined  from  equation  (7-2))  was  81%.  The  sampl- 
ing period  was  set  at  50  gas  meter  resets.  Since  the  gas  meter  was 
calibrated  at  21  mis,  this  entailed  allowing  1.05  liters  of  methane  to 
be  measured   for   each   sample.      The  starting   steady-state    total   VFAs  were 
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Table  7-1.   Algorithm  parameters  for  experimental  runs 


a   =  0.25  parameter  for  calculating  h 
a   =   10  -^  forgetting  factor  parameter 
Pq   =  10'  initial  oovariance  matrix 


"2        = 

0.3 

Hi        = 

0.3 

Hqff  = 

0,1 
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i420  mg  COD/1  of  reactor  volume.   Acetate  and  propionate  were  predominant 
at  concentrations  of  260  mg/1  and  30  mg/1,  respectively. 

Results  of  the  run  are  shown  in  Figures  7-2  to  7-5.  The  total  run 
time  was  10  days.  After  the  initialization  period,  the  algorithm  pre- 
dicted a  minimum  and  increased  temperature.  Temperature  set  point  rose 
from  50°C  to  59°C  while  the  methane  gas  production  increased  32?. 
Consequently,  the  HRT  decreased  from  15  days  to  10  days.  At  this  point, 
the  VFA  concentration  had  only  risen  slightly  to  550  mg  COD/1. 

However,  as  the  temperature  elevated  to  60°C,  the  gas  production 
dropped  dramatically,  which  caused  the  HRT  to  increase  to  30  days.  The 
decrease  in  the  feeding  rate  aided  in  keeping  the  VFAs  from  rising 
appreciably.  At  sample  17,  the  total  VFA  concentration  was  1200  mg 
COD.  Higher  weight  acids  also  appeared  at  this  time. 

Eventually,  the  reactor  restabilized  and  the  gas  production  rose 
and  optimized  at  a  rate  of  approximately  3-3  mls/min  {^.15  liters  of 
methane  per  day).  The  predicted  optimum  was  58°C.  It  should  be  noted, 
however,  that  this  final  rate  was  below  the  starting  value.  This  lower 
production  was  probably  due  to  the  inhibition  caused  by  the  rise  in 
concentration  of  the  higher  weight  VFAs,  such  as  isovalerate  and  iso- 
butyrate.  In  addition,  the  digester  became  very  unstable.  A  one  degree 
increase  in  temperature  done  manually  at  the  end  of  the  experiment  to 
verify  an  optimum  resulted  in  another  collapse  in  the  gas  production. 

Although  the  constant  yield  operating  mode  did  not  prevent  the 
appearance  of  higher  weight  VFAs,  it  did  prevent  the  total  failure  of 
the  digester.  In  fact,  within  three  weeks  the  digester  was  operating  as 
it  had  originally  at  the  start  of  this  run. 
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Figure  7-2.   Experimental  methane  production  for  first  thermophilic  run. 


86 


CO 

■D 


X 


10       15       20       25      30 
Sampling  Interval 


35       40 


Fig-are  7-3.  Experimental  hydraulic  retention  time  for  first  thermo- 
philic run. 
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Actual  reactor  temperature  and  predicted  optimum  tempera- 
ture for  first  thermopnilic  run. 
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Figure    7-5.     Experimental    volatile     fatty    acid    concentration    for    first 
thermopnilic    run.       (Other    acids    are    a   sum   of    isobutryate 
n-butyrate.    isovalerate,    and  n-valerate.) 
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A  probable  reason  for  the  rise  in  VFAs  was  the  effect  of  instan- 
taneous carbon  dioxide  evolution  on  the  gas  meter.  That  is,  the  gas 
meter  measures  methane  accurately  only  when  the  percentage  of  methane 
remains  constant  (see  Appendix  C).  The  acidogenic  population  is  respon- 
sible for  the  majority  of  carbon  dioxide  produced  and  has  substantially 
faster  kinetics  than  methanogens.  Therefore,  a  rise  in  temperature  will 
result  in  a  temporary  increase  in  the  carbon  dioxide  percentage  in  the 
effluent  gas.  This  would  be  especially  unsatisfactory  near  the  methano- 
genic  optimum  where  a  rise  in  temperature  would  cause  a  drop  in  methane 
production  and  a  possible  increase  in  carbon  dioxide  production. 

To  account  for  this  problem  a  second  run  was  started  under  similar 
conditions  with  a  slight  modification.  At  the  beginning  of  each  sampl- 
ing interval,  feeding  rate  was  cut  in  half  for  the  first  five  meter 
resets.  It  was  hoped  that  this  modification  would  slow  the  metabolism 
of  the  acidogenic  bacteria.  This  also  changed  the  reactor  COD  conver- 
sion percentage  to  85it.  As  a  result,  the  starting  total  VFAs  were 
lowered  to  300  mg  COD/1  from  the  above  case. 

The  results  for  this  run  are  shown  in  Figures  7-6  to  7-9.  As  in 
the  above  run,  the  methane  gas  production  rose  dramatically  and  then 
fell  as  the  temperature  rose  to  60°C.  However,  under  the  new  modifica- 
tion of  Initially  cutting  the  feed  rate,  the  VFAs  rose  only  slightly  to 
500  mg  COD/1  with  no  appearance  of  higher  weight  VFAs.  Unfortunately, 
after  sampling  period  21,  a  power  failure  caused  the  reactor  to  shut 
down.  After  restarting  the  reactor,  it  returned  to  normal  operation  in 
24  hours.  It  appeared  that  the  feeding  modification  had  the  advantage 
of  increased  stability.  However,  time  did  not  permit  the  attempt  of  a 
third  run.   The  total  run  time  was  6  days. 
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Figure  7-6.  Experimental  methane  production  for  second  thermophilic 
run. 
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Figure  7-7.  Experimental  hydraulic  retention  time  for  second  th.ermo- 
philio  run. 
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Figure  7-8.  Actual  reactor  temperature  and  predicted  optimum  tempera- 
ture for  second  thermopnilic  run. 
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Figure  7-9.  Experimental  volatile  fatty  acid  concentration  for  second 
thermophilic  run.  (Other  acids  are  a  sum  of  isobutryate, 
n-butyrate,  isovalerate,  and  n-valerate.) 
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7.^.  Mesophillc 

The  digester  was  primed  with  effluent  from  a  gluoose-fed  anaerobic 
digester  operating  at  35°C.  The  yield  was  set  at  9.96  liter  of  methane 
per  liter  of  feed.  The  sampling  interval  was  set  at  50  meter  resets 
(1.05  lit  of  methane).  Using  the  feeding  modification  described  above, 
the  COD  conversion  was  775t.  The  starting  temperature  was  set  at  30°C. 
Starting  VFAs  were  550  mg  COD /I.  Results  of  this  experiment  are  shown 
in  Figures  7-10  to  7-13.   This  run  lasted  approximately  35  days. 

After  the  initialization  period,  the  gas  production  rose  sharply  in 
response  to  a  increase  in  temperature  to  45°C.  Remarkably,  the  methane 
rate  had  tripled  in  the  cime  of  five  days  while  the  total  VFAs  had 
remained  relatively  constant.  Also  of  interest  is  the  decrease  in  HRT 
from  20  to  5  days  HRT.  However,  as  before,  the  temperature  reached  a 
value  (approximately  46°C)  where  the  culture  was  unstable  and  gas  pro- 
duction dropped. 

The  ensuing  elevation  of  VFAs  seemed  to  indicate  the  reactor  had 
failed.  Acetate  rose  from  200  mg/1  to  1200  mg/1  while  propionate  rose 
from  100  mg/1  to  400  mg/1.  In  addition,  isobutyrate  and  isovalerate 
appeared  and  leveled  off  at  150  mg/1  each. 

Nevertheless,  the  algorithm  adjusted  the  temperature  set  point  and 
recovered  the  reactor  within  two  days,  increasing  the  methane  production 
back  to  8.5  rals/min  (12.24  liters  methane  per  day  or  2.0  liters  methane 
per  liter  of  reactor  volume  per  day).  In  contrast  to  the  thermophilic 
runs,  acetate  and  proprionate  decreased  to  300  mg/1  and  120  rag/1, 
respectively.   Isobutyrate  and  Isovalerate  decreased  by  75$. 

At  this  point,  gas  production  began  to  have  an  oscillatory  pat- 
tern.  At  approximately  the  115th  sampling  interval,  however,  the  feed 
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Figure  7-10.  Experimental  methane  production  for  mesophilic  run. 
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Figure  7-11.  Experimental  hydraulic  retention  time  for  mesophilio  run, 
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Figure  7-12.  Actual  rpactor  temperature  and  predicted  optimum  tempera- 
ture for  mesophilio  run. 
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Figure  7-13.  Experimental  volatile  fatty  add  concentration  Tor  meso- 
philio  run.  (Other  acids  are  a  sum  of  isobutryate,  n- 
butyrate,  isovalerate,  and  n-valerate.) 
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tube  clogged  and  was  not  discovered  until  the  119th  interval  (approxi- 
mately 24  hours).  During  this  time,  total  VFAs  fell  to  110  mg  COD/1  and 
the  higher  weight  VFAs  practically  disappeared. 

Once  repaired  and  restarted,  gas  production  continued  its  oscilla- 
tory pattern  imtil  settling  at  a  rate  of  4.1  ml/min  (5.9  liters  methane 
per  day).  The  temperature  converged  to  the  predicted  optimum  of  47.5°C 
at  the  170th  interval.  The  algorithm  maintained  the  reactor  in  this 
state  until  the  l85th  interval  at  which  time  temperature  perturbations 
were  performed  to  verify  the  optimum.  During  the  first  perturbation 
(temperature  increased  to  48.5°C),  the  gas  production  dropped  dramati- 
cally and  higher  weight  VFAs  appeared.  The  reactor  would  not  return  to 
the  operating  conditions  before  the  perturbation  and,  therefore,  a 
perturbation  down  in  temperature  could  not  be  performed. 

7.5.   Discussion 

As  is  the  case  of  many  inaugural  experiments,  more  questions  seem 
to  be  raised  than  answered.  Furthermore,  answers  will  usually  not  match 
theoretical  predictions.  Such  is  the  situation  for  the  experimental 
results  presented  above. 

Foremost,  the  computer  simulations  presented  in  Chapter  3  did  not 
accurately  predict  the  behavior  of  the  anaerobic  digester.  Tnis  is 
invariably  caused  by  the  simple  model  used  to  simulate  the  microbial 
reactor.  In  particular,  inhibition  was  not  considered  in  the  computer 
runs.  Additionally,  the  actual  experimental  implementation  was  not 
accurately  depicted  in  the  simulations.  However,  the  algorithm  did 
succeed  in  increasing  the  production  of  the  digester  although  this 
increase  was  temporary  for  the  thermophilic  runs.   In  the  mesophillc 
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run,  gas  production  reached  extremely  high  levels  before  finally  set- 
tling at  an  optimum.  It  is  apparent  that  much  work  needs  to  be  done  to 
investigate  these  seemingly  transitory  states  and  whether  they  can  be 
maintained. 

Probably  the  most  confusing  part  of  the  experimental  data  is  the 
graphical  representation.  Results  were  plotted  versus  sampling  interval 
number  instead  of  time.  This  was  done  because  of  alternate  sampling 
basis  (Chapter  5).  This  method  replaces  time  with  methane  production  as 
the  independent  variable.  In  fact,  each  sampling  interval  can  be 
thought  of  as  approximately  one  liter  of  methane  produced.  For  example, 
190  liters  of  methane  were  produced  for  the  mesophilic  experiment. 
However,  to  properly  evaluate  the  success  of  the  algorithm,  it  is  neces- 
sary to  know  the  amount  of  time  elapsed.  This  can  easily  be  described 
by  remembering  the  functionality  of  alternate  sampling  basis:  as  gas 
production  increases  the  time  between  sampling  intervals  decreases  and 
vice  versa.  For  the  thermophilic  experiments,  the  time  variation  was 
minimal  with  sampling  occurring  between  three  and  five  times  a  day.  The 
mesopnilic  run,  on  the  other  hand,  exhibited  variations  in  the  sampling 
period  between  two  and  twelve  times  a  day.  Methane  production  for  the 
mesophilic  run  versus  time  is  plotted  versus  time  in  Figure  7-14.  The 
time  variation  of  the  sampling  interval  is  shown  in  Figure  7-15. 

In  addition  to  alternate  sampling  basis,  constant  yield  operation 
(Chapter  6)  causes  some  puzzlement.  In  this  mode,  the  reactor  yield  is 
held  constant  and  HRT  (or  feeding  rate)  is  allowed  to  vary.  This  mode 
has  the  advantage  of  approximately  maintaining  the  environment  of  the 
digester  bacteria  and  isolating  the  effect  of  the  input  variable, 
temperature. 
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Figure  7-1-4.  Experimental  methane  production  for  mesophilic  run  versus 
time.   Compare  to  Figure  7-10. 
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The  results  of  the  thermophilio  runs  were  very  similar.  Both 
seemed  to  show  a  sensitivity  (or  instaoility)  at  or  near  the  optimum. 
Once  passing  this  threshold,  the  gas  production  declined  severely  with 
the  appearance  of  higher  weight  acids.  This  condition  seemed  to  be 
irreversible  as  the  algorithm  was  not  able  to  recover  the  digester  to 
the  original  state.  Nevertheless,  the  constant  yield  operating  mode 
appeared  to  save  the  digester  from  total  failure  keeping  it  in  an  oper- 
able condition.  The  alteration  of  delaying  the  feeding  after  each 
temperature  change  appeared  to  help,  although  the  effect  could  be  due  to 
the  increased  yield  (or  conversion  percentage).  Decreasing  the  head 
space  of  the  reactor  would  .'^educe  the  delay  involved  in  implementing 
constant  yield  policy  (Appendix  C).  This  would  most  likely  improve 
operation  of  the  digester  at  these  thresholds. 

Instability  at  high  methane  production  was  similarly  found  in  the 
mesophilic  run.  However,  the  digester  behavior  was  distinctly  different 
from  thermophilic.  The  sharp  decrease  in  gas  production  and  subsequent 
rise  in  VFAs  was  temporary  and  the  reactor  was  easily  recovered.  Fur- 
thermore, the  gas  production  seemed  to  oscillate  before  settling  at 
U7.5°C.  This  phenomenon  may  be  due  to  the  adaption  of  the  bacteria  to  a 
new  metabolism  that  is  better  suited  for  the  new  temperatures.  Another 
reason  could  be  population  shifting. 

Of  particular  interest,  is  the  range  of  operation  for  the  meso- 
philic run.  In  fact,  it  is  debatable  whether  the  digester  was  meso- 
philic or  thermophilic  since  the  operating  range  for  the  last  half  of 
the  r'un  was  in  the  region  generally  accepted  to  be  intermediary.  This 
result  supports  the  assertions  of  Chen  et  al.  [59]  that  a  smooth  transi- 
tion exists  between  MO  and  50°C  (Chapter  4). 
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In  conclusion,  the  experimental  results  were  successful  in  two 
areas.  Firstly,  the  mesopnilic  digester  methane  production  was  in- 
creased by  75  per  cent  over  the  initial  starting  state.  This  supports 
the  further  experimentation  of  adaptive  optimization  on  microbial  growth 
systems.  Secondly,  the  experiments  revealed  many  areas  of  future  re- 
search. Among  these  are  the  investigation  of  digester  behavior  at 
optimal  conditions  and  the  use  of  adaptive  algorithms  to  recover  inhi- 
bited digesters. 


CHAPTER  8 
SUMMARY 

Adaptive  optimization  was  presented  in  this  work  as  a  method  of 
maximizing  production.  The  advantages  offered  by  this  approach  over 
other  methods  are  many.  Primarily,  continuous  reactors  can  be  optimized 
quickly  with  little  or  no  advance  knowledge  of  the  process.  Further- 
more, optimization  is  performed  on-line  and  thus  avoids  the  problems 
associated  with  ad  hoc  experiments.  However,  a  substantial  effort  is 
needed  to  properly  tune  these  algorithms  and  realize  their  benefits. 

It  is  apparent  that  much  of  the  problems  with  this  approach  stem 
from  measurement  error.  These  errors  affect  the  parameter  convergence 
and  subsequently  impair  the  ability  of  the  algorithm  to  successfully 
predict  an  optimum  input.  Better  results  would  be  obtained  if  the 
adaptive  optimization  algorithm  were  designed  to  handle  such  problems. 
That  is,  numerical  simulations  should  be  done  with  actual  process  condi- 
tions in  mind  (e.g.,  measurement  error,  controller  resolution,  and 
upsets) . 

It  would  also  be  of  some  good  consequence  to  perform,  if  possible, 
simple  sensitivity  studies  to  better  understand  the  limitations  of  the 
controller  input  (such  as  maximum  and  minimum  input  change  limits). 
These  experiments  would  necessarily  breach  the  advantage  of  no  a  priori 
knowledge,  but  could  play  a  big  part  in  the  success  of  the  algorithm. 
However,  these  experiments  could  easily  be  designed  to  be  included  in 
the  initialization  period. 
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The  innovations  of  constant  reactor  yield  operation  and  alternate 
sampling  basis  revealed  some  interesting  results.  Unusually  high  gas 
production  was  achieved  in  the  mesophilic  run  with  relatively  low  vola- 
tile acid  concentration.  This  was  mainly  due  to  the  constant  yield 
operation  mode.  This  mode  offers  some  promise  in  improving  the  optimi- 
zation of  other  continuous  microbial  growth  processes.  In  addition,  the 
method  of  alternate  sampling  basis  showed  success  in  reducing  the  time 
between  samples  and  thus  improved  the  performance  of  the  adaptive  opti- 
mizer. 

Experimental  ri.ms  on  an  anaerobic  digester  exhibited  some  unique 
characteristics  of  a  mixed  culture  system.  Instability  at  the  optimum 
seemed  to  be  present  in  both  thermophilic  and  mesophilic  reactors.  This 
indicates  that  it  could  be  better  to  operate  the  digester  just  below  the 
optimum. 

Finally,  much  Information  was  revealed  by  the  raesophilii^  run.  The 
oscillations  could  be  the  result,  in  part,  of  the  shifting  of  bacterial 
populations.  Information  in  population  dynamics  could  be  obtained  by 
repeating  these  experiments  and  enlisting  the  use  of  biosensors. 


APPENDIX  A 
JUSTIFICATION  OF  STEP  SIZE 

From  step  7  of  the  algorithm  of  section  3.2  we  have 


D(t)  =  h(t)  D(t-1)  +  (1-h(t))  D*(e(t))  +  V(t)        (A-1) 

Since  the  estimates  6(t)  are  not  perfect  there  is  some  error 
D(t)  in  D»(e(t)) 

D(t)  =  D»(d(t))  -  D»(t))  (A-2) 

where  9(t)  are  the  correct  (unknown)  parameters.  In  view  of  equation 
(A-2) 

D(t)    =   h(t)    D(t-l)    -    (l-h(t))    5(t)    +    (l-h(t))    D»(e(t))    +   V(t)      (A-3) 

The  first  terra  of  the  right  hand  side  can  be  considered  as  a  caution 
term  while  the  second  term  is  the  error  in  D(t).  It  is  reasonable  to 
limit  this  erroneous  part  to  a  small  fraction  6  of  the  caution  part, 
i.e. 

(l-h(t))  I  5(t)  I  =  5  h(t)  D(t-1  ) 


or 
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|D(t)|  +  6  D(t-1 ) 

To  obtain  an  expression  for  h(t)  in  terms  of  known  quantities  an  esti- 
mate for  |D(t)|  is  needed. 


Since 


D»(e(t) 


c(t)  -  qs  a(t)  fA-S') 


where 


a(t)  =  1  -  -Zq  a.(t) 


b(t)  =  .Zq  b.(t) 


9D*(e) 

D*(8(t))  -  D*(8(t))  +  —TK—    I-    ^9(t)  -  9(t) 
-  -  dQ  e(t) 


which  gives  that 


wh 


q3°a(t)  -  2D*(9(t))b(t)  -  o(t) 

5(t)  =  X 'A-6) 

2b(t) 


ere  i(t)  =  a(t)  -  a(t),  b(t)  =  b(t)  -  b(t),  o(t)  =  c(t)  -  c(t) 


After  the  parameter  estimates  have  been  updated  the  a  posteriori  estima- 
tion error  v(t)  can  be  calculated  from  equation  (3-6).  Since  typically 
biochemical  processes  have  slow  dynamics 
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x(t-1 )  =  x(t-l-i),  i  =  l,...,n  and  D(t-1  )  =  D(t-l-i),  i  =  1 m, 


which  in  view  of  equation  (3-2)  gives 


v(t)  ~  a(t)  x(t-1)  *   b(t)  D(t-1)  +  c(t)  (A-7) 


This  suggests  that  very  rough,  order  of  magnitude  estimates  are 


a(t)  -^^'^ 


x(t-1  ) 


5(t):^{^  (A-8) 


G(t)  -  v(t) 


where  v  is  as  in  equation  (3-8),  i.e.  we  filter  v  whenever  it  decreases 
in  order  to  protect  against  a  low  value  due  to  error  cancellation. 
Substituting  (A-8)  in  (A-6)  and  the  outcome  in  (A-U)  we  obtain  equation 
(3-7)  with 

a   =   26D(t-1)|b(t)|       ^  (^_g) 

o        2D*(e(t)) 

qs 
— a +  1  +  


x(t-1)         D(t-l) 


APPENDIX    B 
THE   PARAMETER    ESTIMATOR 


Defining 


8^=    Ca^,    a^,....    a^.    b^ , . . . .    b^,    o,    b^] 


i   (t-1)    =    [x(t-1) x(t-n),    D(t-2),...,    D(t-m),    1,    D(t-1)] 


model  (2)  is  rewritten  as 


x(t)  =  /(t-1  )  6  (3-1  ) 


The  estimator  is  then  as  follows 


Initialization 


1(0)  =  9        (If  no  a  priori  estimate  is  available  the 

'°  T 

typical  choice  is   9  =  [0 0,1]) 


U(0)  =   [u  fO),  u^(0) u  (0)]  =  I 

-1     -2         -n 


^i(O)  =  1/Y  ;  i  =  1 n 
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Ill 


Recursive   steps 


f   =   U'^Ct-l  )    ti;(t-1  ) 


rj=  L^y  r^ f^] 


v^=   dj (t-1  )f^;    i  =  1  ,  ..  .,n 


w   =    1    +      Z   V   f 
1  =  11    i 


T 
e   =   x(t)    -  ^   (t-1)    6(t-l)        (a  priori  error) 


2 

A   =   1    -  £   /[(l+w)o]  (forgetting   factor) 


If  X   <  0,7   then   A   =   0.001  (covariance  resetting) 


'r  '  *  ^1^1 


6j:=   B^_^+  v^f^;    i    =  2 n 


d^(t)    =  d^  (t-1  )/(B^A) 


d^(t)    =   dj,(t-l)    6^_^/(B^A);    i    =   2,...,n 


«.^=  -  ^i/B^_^  ;   i    =2 n 


k2  =    [v^  ,    0,    ...,    0] 
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k,    ,=   k   +   v^    u. (t-1  );    i    =   2 n 

-i+1      -1        i    -1 


-n+1      n 


8(t)   =   e(t-1  )   +   k  e  (new  estimate) 


u^(t)    =   u^(t-l) 


u    (t)   =   u.(t-l  )    +   i^k^;    i   =   2,...,n 


APPENDIX    C 
MATERIALS  AND   METHODS 

C.I .     Experimental   Set-up 


Reactor  Construction 

The  bioreactor  used  for  the  digestion  process  is  constructed  from 
6"  ID  wrought  iron  pipe  coated  with  a  water  base  epoxy.  Two  1"  baffles 
are  attached  to  the  walls  to  assist  circulation  of  biomass.  A  stirring 
rod  with  three  propellers,  also  coated  with  epoxy,  is  operated  continu- 
ously at  150  rpm  (Figure  C-1).  The  total  working  volume  of  the  reactor 
is   ~  6.5   liters.      The  liquid  volume   is  6.2   liters. 

Gas  Measurement   System 

The  gas  measurement  system  is  a  raonoraeter  type  apparatus  developed 
in  the  Bioprocess  Engineering  Research  Laboratory  at  the  University  of 
Florida. 

The  device  is  constructed  from  1  .25"  poly  vinyl  chloride  pipe  and 
fittings  into  a  U-shape  (see  Figure  C-2).  The  monometer  fluid  is  stand- 
ard IOW-3O  motor  oil.  A  float-switch  level  sensor  is  attached  to  the 
end  of  the  U-tube  open  to  the  atmosphere.  A  two-way  valve  is  attached 
to  the  other  end,  normally  open  to  the  gas  line  of  the  reactor  and 
normally  closed   to   the  ambient  air. 

Typical    operation    of    the    gas    meter    is    as    follows.       When    the    gas 

pressure    in    the    reactor    becomes    sufficient    to   displace   a    predetermined 

volume  of  oil   in   the  gas  meter,    the  float  switch  activates  a  relay 
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Figure  C-1.   Schematic  of  experimental  anaerobic  digester 
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timer.  For  the  following  three  seconds,  the  effluent  gas  line  is  closed 
and  the  gas  side  of  the  meter  is  opened  to  the  atmosphere,  allowing  the 
manometer  to  return  to  equilibrium. 

The  meter  is  calibrated  by  injecting  nitrogen  gas  into  the  effluent 
gas  line  until  the  float  switch  trips  twice  and  measi-iring  the  volume 
injected  between  the  first  and  second  trips.  This  trip  volume  can  be 
reduced  or  increased  by  adding  or  removing,  respectively,  monometer 
fluid. 

The  gas  meter  is  connected  to  the  reactor  with  tygon  tubing  (Cole 
Palmer).  An  adsorber  filled  with  indicating  soda  lime  is  placed  in- 
line, close  to  the  reactor,  to  remove  carbon  dioxide,  so  that  methane 
production  can  be  measured.  A  metal  line  is  used  to  exit  the  reactor  to 
bring  the  gas  to  ambient  temperature.  A  quantitative  description  can 
also  be  derived  to  analyze  measurement  errors: 

The  gas  system  is  divided  into  three  volume  sections  (see  Figure  C- 

2).   Section  three  is  further  divided:   V-,  =  V,  +  V„  and  V„  =  V„,  +  Vn 

3    1-1    m       m    ma    u 

where  V^  is  the  line  volume,  V  is  the  meter  volume  (volume  after  two 
way  valve),  V^^^^  is  the  meter  volume  at  atmospheric  pressure,  Vn  is  the 
oil  displacement  volume.  As  described  previously,  when  Vn  reached  V  , 
critical  displacement  volume,  the  float  switch  will  activate  and  the 
production  of  the  calibrated  trip  volume,  V-p,  will  be  registered. 
Production  rate  is  measured  by  dividing  V~  by  the  time  between  trips. 
For  Vj^Q  >  Vq  >  0,  the  relationship  between  system  pressure  and  volume  is 


P  g 

P  =  P  +  ^  V^  (C-1) 

3    a    A   D 
m 
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Figure  C-2.      Schematic  of  gas  meter, 
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where  P^^  is  the  density  of  the  oil  and  A^  is  the  area  of  the  meter  tube 
cross  section.  Assuming  ideal  gas  law,  the  state  equations  for  each 
section   is: 


n.  =  (P   *  BV^)V,/RT,  (C-2) 

1     a     D  1    I 


n_  =  (P  +  ev^)V,/RT  (C-3) 

2     a     D   2   a 


n,  =  (P   +  SV^)(V,  ^  V   +  V^)/RT       (C-4) 
3     a     D   L    ma    D    a 

ZP  g 

Here,  3  is   "  ,  n  number  of  moles  is  section  i,  P^  is  ambient 
m 
pressure,  T  is  ambient  temperature,  T,  is  the  temperature  of  the  reac- 
tor. 

The  gas  flow  is  assumed  to  be  positive,  i.e.,  from  reactor  to  gas 
meter,  and  to  undergo  an  immediate  temperature  change  to  ambient  from 
section  1  to  section  2.  (Note  that  negative  gas  flow  might  be  possible 
because  of  pressure  and  temperature  changes.  These  special  cases  will 
be  handled  later.)  The  assumption  will  be  made  that  ambient  temperature 
and  pressure  are  constant.   The  mass  balances  for  each  section  is 


dn^ 

^=^G-"l2  ^^'^^ 

dn^ 


dn. 
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Here  r^  is  gas  production  rate,  n. .  is  molar  flow  rate  from  section  i  to 
section  j  ,  x  is  mole  fraction  of  methane  in  section  2.  It  should  be 
noted  that  the  complete  adsorption  of  methane  from  section  2  to  section 
3  is  a  very  accurate  description.  In  trial  runs,  carbon  dioxide  adsorp- 
tion was  measured  at  99.9$  for  gas  rates  of  50  rals/min  of  60$  carbon 
dioxide  and  H0%  methane.  Typical  flow  rates  for  the  anaerobic  digester 
used  were   4   to  20  mls/min. 

Combining  equations   C-5   to  C-7  gives: 


dn         dn.  dn, 

— -  +  — -  +  -  — ^  =  r  (C-8) 

dt  dt        x     dt  G 


A  methane  balance  gives  the  equation  for  mole  fractions  in  sections 
1  and  2  (assuming  complete  mixing  in  both  sections): 


dt     T„     2         dn 

^  =  i^iy).       =!i  (c-10) 

dt     T,  I    r^ 

I  Li 


Where  y  is  the  mole  fraction  in  section  2  and  z  is  the  mole  fraction  of 
gas  generated,  Vq. 

By  inspection  of  equations  C-9  and  C-10,  it  should  be  noted  that  in 
order  to  reduce  response  time  due  to  a  change  in  methane  production 
rate,  zrQ,  the  time  constants,  x, ,  of  sections  1  and  2  should  be  made  as 
small  as  possible.  This  would  be  facilitated  by  reducing  the  volumes  V, 
and  V2.  However,  V2  cannot  be  reduced  to  the  point  that  the  adsorption 
of  carbon  dioxide  affects  the  concentration  of  carbon  dioxide  in  the 
reactor  head  space.   A  limiting  distance  between  the  adsorber  and  the 
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reactor    can    easily    be    approximated    from    molecular    diffusion    laws    as 
shown  by  Treybal   [68]: 

X,    .      =  ^5_  ^n   [1^]  (C-11) 

'"^^        V    . 
ram 

where   A   is    the  area  of    tubing   in  section   2,    D   is  diffusivity  of  methane 

and    carbon    dioxide,    V    .      is    the    minimum    volumetric    flow    rate    of    total 

min 

» 
gas,    X   is    the  mole    fraction   of  methane  at  S,    .    ,    x     is    the  mole   ^paction 
°  min 

of  methane  gas   production. 

At  25°C  and  1  atm,  the  di*^fusivity  is  approximately  1  ,000  mls/m-min 
[65].  For  a  minimum  flow  rate  of  1  ml/min  and  a  change  in  x  equal  to 
0.5?  of  X    ,    the  minimum  distance  for    the  system  used  here   is    ~    10  cm. 

Assuming  the  system  concentrations  are  in  equilibrium,  i.e.,  x=y=z, 
equations  C-2  to  C-4  can  be  differentiated  and  substituted  into  equation 
8: 

P   V      +    3V.V_      dT, 
(•    a   1             1    Di        1 
r,    +   ' 


G        ^  2  ^      dt   " 

R  lyy  "  Tri ^  "dt  ^^  '/^ 

At  most  times  —-rr   will  be  zero  except  when  the  reactor  temperature 
changes,  in  which  case  it  will  cause  positive  errors  in  gas  production 


dT 


measured  for  — --  >  0.   An  estimation  of  this  error  can  be  made  by  assum- 


dT. 


ing  the  term  in  front  of  — —  in  equation  12  is  constant.  This  reveals 
that  a  1  °C  change  up  in  temperature  simulated  the  production  of  3-35  x 
lO"-*  moles  of  gas  or  1  ml  of  methane.  For  the  sampling  period  used  in 
this  experiment  this  amounts  to  a  0.1  J  error. 
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If  no  temperature  change  is  assumed,  then  equation  12  can  be 
integrated  to  get: 

2 

P  t  =  ^[-i  -  — - ^ ^  V  *  ^]  -  C        (C-13) 

"^G^   R^T,  yT  D   T  ^ 

I  a  a 

The  volume  needed  to  reset  the  gas  meter  is  V^^.  When  this  volume 
is  reached,  the  meter  is  closed  to  the  effluent  gas  line  and  opened  to 
the  atmosphere.   The  moles  released  from  the  system  is: 


An  = 


(P   +  BV   )V   +  6V^ 

^  a   '^  mo^  DC   "  Dc  (C-14) 


RT 
a 


At  steady-state,  a  molar  balance  gives  An  =  i"  At,  where  At  is  the 
time  between  gas  meter  trips,  or  for  equation  C-13  the  time  for  Vp  to  go 
from  Vp  ,  initial  volume,  to  V^^.  Note  that  the  gas  measurement  will  be 
highly  sensitive  to  changes  in  ambient  temperature  and  pressure.  Also, 
Vq  is  greater  than  zero  because  of  the  high  pressure  (P  =  P^  +  ^^dc^  ^" 
the  reactor  and  gas  lines  when  the  gas  meter  goes  back  to  the  normally 
open  position.  V^^  can  be  found  by  substituting  C-14  into  C-13  or  by 
doing  a  molar  balance  of  the  system  at  the  point  when  the  gas  meter  is 
reset. 

Anaerobic  Overflow 

An  anaerobic  overflow  is  used  to  maintain  the  liquid  level  and 
anoxic  conditions  in  the  overflow  line  (Figure  C-1).  The  effluent  line 
is  drawn  from  the  three  liter  mark.  A  line  from  the  head  space  of  the 
reactor  is  attached  to  "y"  connector  kept  at  the  liquid  level  of  the 
system.   One  branch  of  the  connector  is  attached  to  the  effluent  line 
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and  the  other  branch  is  attached  to  the  overflow  line.  The  end  of  the 
overflow  line  is  submerged  in  a  side  arm  flask  filled  with  distilled 
water.  In  this  apparatus,  excess  fluid  in  the  reactor  overflows  at  the 
"y"  and  flows  into  the  side  arm  flask.  Fluid  from  the  sidearra  is  col- 
lected in  a  plastic  jug  and  discarded.  A  pimp  is  used  to  recirculate 
fluid  from  the  effluent  line  back  into  the  reactor.  This  is  done  in 
order  to  keep  biomass  from  settling  and  inhibiting  flow  through  the 
line. 

Temperature  Control 

The  temperature  of  the  system  is  controlled  via  an  Omega  CN200 
controller.  A  platinum  RTD  is  used  to  measure  temperature  with  an 
accuracy  of  0.1°C.  Heating  is  achieved  through  two  50  U  heating  slabs 
connected  in  series  on  opposite  sides  of  the  outside  of  reactor.  Heat 
input  is  adjusted  by  changing  the  percentage  time  full  power  is  applied 
to  the  slabs.  Cooling  is  achieved  by  an  internal  coil  constructed  of 
stainless  steel  tubing.  A  water/ethylene  glycol  mixture  kept  at  0°C  is 
circulated  through  the  tubing  to  lower  the  temperature.  Fiberglass 
insulation  is  used  to  reduce  heat  flow  to  the  ambient. 

Feed  System 

The  feed  is  pumped  through  a  valveless  micro  piston  metering  piamp 
(Fluid  Metering  Incorporated).  Both  piston  and  cylinder  housing  are 
constructed  of  ceramic  The  cylinder  is  0.25"  ID  and  the  piston  is 
ground  at  the  factory  for  airtight  fit.  Two  rulon  seals  are  placed  on 
the  piston  shaft  to  insure  fit. 
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The  feed  rate  can  be  adjusted  in  two  ways.  The  stroke  length  of 
the  piston  can  be  adjusted  on  the  pump  head  to  alter  the  amount  of  fluid 
delivered  per  revolution.  Feed  rate  can  also  be  adjusted  by  changing 
the  speed  of  the  driver.  The  driver  is  a  variable  speed  1/10  HP  DC 
motor  (Cole  Palmer)  and  is  adjusted  by  a  solid  state  speed  controller. 
Rotation  speed  can  be  monitored  by  a  computer  through  a  cam  attached  to 
the  motor  drive  and  a  rocker  arm  relay  switch  (see  Electronics). 

The  pump  head  is  calibrated  before  each  experiment  to  determine 
volume  per  revolution.  The  calibration  is  checked  throughout  the  exper- 
iment by  measuring  the  volume  pumped  from  the  feedtank. 

An  0.065"  ID  TFE,  0.125"  00  teflon  tubing  is  used  to  deliver  the 
feed  from  the  tank  to  the  pump  and  from  the  pump  to  a  top  port  on  the 
reactor.  The  tubing  is  coupled  to  the  pi-imphead  with  1/8"  teflon  com- 
pression fittings.  The  tubing  is  also  connected  to  the  reactor  by  a 
1/8"  stainless  steel  compression  fitting. 

A  70  watt  UV  light  is  used  to  keep  the  feed  lines  free  of  cell 
growth.  This  was  necessary  to  keep  the  tubing  and  p'omphead  free  from 
clogging.  Tubing  is  cleaned  with  hot  water  and  ethyl  alcohol  before 
each  experiment. 

Computer  Implementation 

An  IBM  compatible  XT  personal  computer  is  used  to  control  and 
monitor  the  overall  experiment. 

Communication  with  the  temperature  controller  (monitoring  tempera- 
ture, adjusting  set  point)  is  achieved  via  an  RS232C  port. 

Gas  measurement  is  calculated  by  dividing  the  calibrated  meter 
volume  with  the  time  between  float  trips,  monitored  by  the  computer  via 
a  digital  input  line  and  signal  debouncer. 
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The  piomp  is  turned  on  and  off  through  a  digital  output  line  and  a 
110V/5V  relay.  The  computer  monitors  the  rotations  through  a  cam  relay 
and  digital  input  line. 

C.2.   Feed  Stock  and  Preparation 

The  digester  is  fed  continuously  on  a  sterile  feed  o^"  glucose  based 
on  a  defined  medium  reported  by  Owen  et  al.  [66]  supplemented  with  tung- 
stic  acid  and  nickelous  chloride  (Table  C-1). 

The  feeds  are  prepared  by  dissolving  glucose,  yeast  extract  (Dif- 
co)  ,  cassamino  acids  (Difco),  and  sodium  bicarbonate  in  I960  mis  of 
distilled  and  deionized  water  in  a  2  liter  polypropylene  flask.  Nitro- 
gen is  bubbled  through  the  mixture  for  10  minutes.  The  '"lask  is  then 
closed  with  a  silicon  rubber  stopper  fitted  with  three  glass  tubes,  each 
fitted  with  silicon  tubing.  The  stopper  is  fastened  to  the  flask  with  a 
hose  clamp  and  wire  to  permit  autoclaving  without  losing  fluid.  Each 
tube  is  covered  with  aluminum  foil  wrap  and  closed  with  tube  clamps. 

The  feed  is  then  autoclaved  for  20  minutes  at  260°F.  The  feed  is 
removed  and  allowed  to  cool  in  an  ice  bath. 

Tnen  30  mis  of  solution  S4  and  10  mis  of  solution  S3  are  added  to 
the  feed  via  a  sterile  syringe. 

The  feed  line,  a  0.125"  CD  teflon  tube,  is  inserted  into  one  of  the 
glass  tubes.  A  sterile  nitrogen  bag  is  supplied  to  an  'inused  glass  tube 
for  suction  relief  caused  by  pumping  feed  to  the  reactor.  The  feed  is 
kept  at  3°C. 
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C.3'      Analysis 
Volatile   Fatty  Acids 

The  volatile  fatty  acid  concentrations  were  obtained  by  gas-liquid 
chromatography.  Samples  were  acidified  with  20?  (V/V)  H^PO^  to  a  final 
concentration  of  2%  (V/V).  The  samples  were  then  centrifuged  in  poly- 
propylene tubes  to  remove  precipitate.  Tr\e  supernatent  was  injected 
onto  1.8  m  X  2  mm  glass  column  of  ^0%  SP-1000  on  100/120  chromosorb 
W/AW  (Supelco,  Inc.).  The  column  was  maintained  at  140°C  with  nitrogen 
as  the  carrier  gas  at  40  mls/min.  The  injector  was  set  at  160°C.  The 
flame  ionization  detector  was  set  at  200°C.  Volatile  fatty  acids  ana- 
lyzed were  acetic,  propionic,  isobutyric,  n-butyric,  iso-valeric,  and  n- 
valer ic. 

Gas  Composition 

Gas  composition  was  determined  by  gas  chromatography.  The  gas 
chromatograph  was  equipped  with  two  stainless  steel  columns  operated  at 
55°C  and  a  thermal  conductivity  detector.  The  <"irst  column  was  2  m  x 
3.2  mm  packed  with  80-100  mesh  Porapak  Q  (Supelco,  Inc.),  and  the  second 
column  was  3.35  m  x  4.8  mm  packed  with  6-80  mesh  molecular  sieve  13X 
(Supelco,      Inc.).  Helium     was      the     carrier     gas,      maintained     at     30 

mls/min.        Gases    monitored    were    nitrogen,     methane,     carbon    dioxide    and 
oxygen. 

COD   Analysis 

Chemical  oxygen  demand  of  yeast  extract  and  cassamino  acids  was 
determined  by  closed  ref lux-titrimetr ic  standard  method  [67]  (HACH 
Reactor    Digestion    Method,     EPA   approved).       The   conversion    factors   were 
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determined  to  be  0.92  g  COD/g  cassamino  acids  and  1.71  g  COD/g  yeast 
extract.  Using  these  values  and  the  theoretical  value  of  1.067  g  COD/g 
glucose,  the  feed  COD  level  was  calculated  to  be  38.89  g  COD/lit  feed. 
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Table  C-1.   Medium  ^or   glucose  digester 


Component 


Stock  Solution 

mg/1 


Final  Concentration 
mg/1 


Glucose 
Yeast  Extract 
Cassamino  Acids 
NaHCOo 
S3 


(NHi^)2HP0n 


CaCl2-2H20 
NH4CI 


MgCl^-eH^O 

KCl 


MnCl  -^H  0 
CoCl  -BH  0 
H3BO3 

CuCl^-ZH^O 
Na2Mo0j^'2H   0 
ZnClo 


NiCl    -SH  0 


H2W0i^ 


14,420 

15,030 

23,940 

10,800 

78,030 

1  ,197 

1  ,800 

342 

162 

153 

126 

135 

6.3 


28,800 
3,200 
3,200 
6,720 

72.1 

225 
359 

1  ,620 

1,170 
18.0 
27.0 
5.13 
2.43 
2.30 
1.89 
2.0 
0.1 
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